{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The objectives of the lab\n",
    "\n",
    "The purpose of this lab is to reproduce tables 3.1 and 3.2 from the third chapter of the book \"Elements of Statistical Learning\" from Hastie, Tibshirani and Friedman, as shown bellow. The objective of the homework is to reproduce table 3.3.\n",
    "\n",
    "## 1. Prepare the data : \n",
    " \n",
    " Raw data is available on line, download it from  moodle ({\\tt Data.txt} file) or from the web at \n",
    "https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "url_data = \"https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data\"\n",
    "df = pd.read_csv(url_data, delimiter='\\t') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For some reason, data has to be normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "variables = df.columns[1:9]\n",
    "df[variables] = df[variables].apply(lambda x: (x - x.mean()) / x.std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split the data into training and test sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training set : n = 67 samples and p = 8 dimensions\n",
      "Test set : n = 30 samples and p = 8 dimensions\n"
     ]
    }
   ],
   "source": [
    "# Get the training and test sets\n",
    "Y_train = df.loc[df[\"train\"]==\"T\", 'lpsa'].to_numpy()\n",
    "X_train = df.loc[df[\"train\"]==\"T\", variables].to_numpy()\n",
    "print(\"Training set : n = {} samples and p = {} dimensions\".format(X_train.shape[0], X_train.shape[1]))\n",
    "\n",
    "Y_test = df.loc[df[\"train\"]==\"F\", 'lpsa'].to_numpy()\n",
    "X_test = df.loc[df[\"train\"]==\"F\", variables].to_numpy()\n",
    "print(\"Test set : n = {} samples and p = {} dimensions\".format(X_test.shape[0], X_test.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Compute the correlation of the predictors in the cancer data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fd97031c5d0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD8CAYAAABTjp5OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd3wURf/H37NX00lyaRBAepMmHVGkKeCjyKNi5VEURBBFmtj1EcX2iIgNKfaKHRVE6V2KUqUjNb23S+5ud35/bMzlktwlFEXy2/frda/kdr4z89nZ2e/Nzk4RUkoMDAwMDP75KOdagIGBgYFBzTActoGBgcF5guGwDQwMDM4TDIdtYGBgcJ5gOGwDAwOD8wTDYRsYGBicJxgO28DAwMAPQoi3hRBpQohdfsKFEGKWEOKgEGKHEOKicmEDhRD7SsMePBt6DIdtYGBg4J93gYEBwgcBzUo/dwFvAgghTMDrpeGtgZuEEK3PVIzhsA0MDAz8IKVcDWQFMBkCvC91NgJ1hBAJQFfgoJTysJTSBXxaantGmM80gepwZxw+51MpC++981xLAEB6znlRkLIt+FxLAOCCqS3PtQSIiT/XCgBQ12881xIQYf+MehHy6IfiTNM4FZ9jjWkyGr1l/CdzpJRzTiG7esDxct9PlB6r6ni3U0i3Sv5yh21gYGDwT6XUOZ+Kg65IVT8wMsDxM8Jw2AYGBrULTf07czsB1C/3PRFIAqx+jp8RRh+2gYFB7UL11Pxz5iwE/lM6WqQ7kCulTAY2A82EEI2EEFbgxlLbM8JoYRsYGNQqpNTOWlpCiE+AywCHEOIE8ARg0fORs4FFwGDgIFAEjCgN8wghxgFLABPwtpRy95nqMRy2gYFB7UI7ew5bSnlTNeESuMdP2CJ0h37WMBy2gYFB7eIstrD/aRgO28DAoHbx9750/FsxHLaBgUHtwmhhGxgYGJwfyLMz+uMfyd/usB+dPoPV6zYRFVmHbz6cXSlcSsmzM2ezZsNm7HYbzzwyidYtmgKwduMWnps5G1XTuPaqgYwcPgyA3Lx8Jj32LEkpqdSNj+OlaQ8RER4WUIe5fReC/jMOFBOuFT9QsvAT3/BOFxM0bARoEqmpON9/DXXfroBxRUgYweMfR3HEo2WkUPTKf5GFBf41dOhK8Ag9nZJlP1Dyzcc+4dZe/bFdU/rOo9hJ0dyXUY8eAsB25XXY+l0JEtRjhyl843lwu/SwgUOxDRoKqor71404P3wrYFmEXNKJuEdHI0wKOQuWkDnn80o2cY+NJrR3FzRnCclTZ1D8u65DCQshYfp4bM0aApLkB2fi3LaXmPuHE9qvO0gNT2YuyVNn4EnzP8N33R9pvLDsdzQpGdquPnd0a+oTvvlYJhO+3kLdCH1GXr/m8Yzu2QyAj7b+wVc7jiEl/LtdA27t3AiA19fuY+WBVIQQRAVbeWpwe2JD7QHLYt2+E7zw7UY0qTG0awvu6NO+ks3mQ8m8uHAjHk0jMtjO/DFX6jrW7uKrX/YhgX93bcGtl1wIwIzvN7F6zzEsJoXE6HD+O+wSwoNsfjWYmnfEdvUdIBTcm5fiXvl1lXZKYlOC7nmW4o9noO7cgIiIxnbDfShhkUip4fnlZ9zrfiizt/QcjKXnIKSmou7ZimvxBwHLwtS4HdYrhoNQ8GxbiXv9d1XrSGiMfcSTlHz1KurezQCYuw7E0vEykBIt/QQlC+eA6sbS+zrMzS9CSglFeZQsfAtZkBNQxylzFl86/tP42x32NYMHcPO1V/PwtP9VGb5mw2aOnUhi0Wfz2bF7L9P+9xqfzJ2Jqqo8/dLrzJ05nfhYBzeMHE+fXt1o0qgh8z5YQPfOHRg5fBjzPljA/A8XMHFsgOnoQiFoxHgKp09By0wn7JnZuLeuRzt5tMzEs2sr+VvXAaA0aEzIfU+QP/m2gHFtQ27Gs+tXShZ+gu3qm7BdfTPFn/iZRKUoBN85noJpk9Gy0gl7djbuLevQTng1qGnJFDwxHllYoDv30ZPIf3gsIsqBbfC15E24DVwuQiY8gfXivrhW/oi5TQcsXXqRN+lO8LgR4XUCXxBFIf7JsRy7/RHcKRk0+nIm+cs34jronVUb0rsz1ob1ONR/JPYOLYh/ahxHrpsAQNyjoylcvZWT904HixnFrjuizHlfkD5TdwiR/7kax7ibSXn8tSolqJrk2Z93M3tYN+LC7NzywVp6N4mjicP3R7djYhSvXtvF59jB9Hy+2nGMD2/thcUkuOfzTVzSJJaGkSHc1qUx9/RqAcDHW/9gzvoDPHp5W79FoWoaz369ntmjBhIXEcItry6kd+sGNImLLLPJc5bw7Nfref3OK0iIDCWrwKnrSMniq1/28eG9Q7CYFO6Zv4RLWtanYUwE3ZvX5b5BnTGbFGYu2sTbK7Zz/+CuVYsQCrZrRuGc919kbiZB417A8/tmZNqJSnbWQcNR92/zHtM0XN+/h5Z0GKx2gu/7H54D25FpJzA1vhBT6y4UvTwBVA8iJMJvOejpC6yDbqP4o+eQeVnY73wKz/6tyIykynb9bkA9vMN7KCwSS9fLcc6eCh43tn/fi7lNdzw71uDe8APuVV8AYO5yOZZLhuJa/E5gLadKLe4S+dsnznTu0DZg63fF2o1cPbAfQgjaX9iK/PwC0jOy2LlnPw0S61K/XgIWi4VB/XqzfI2+BsOKNRsYMqg/AEMG9Wf56g0BNZiatkRLSUJLSwbVg2vDciydL/Y1Kiku+1fY7Pw5qzRQXEunnrhWLwHAtXpJ5TQraTipp+Px4F63HGsFe3X/7rIWunrgd5ToGK8mxYSw2kAxgc2OlpUBgO3yIRR/8zF43ADIvMCtl6B2zXEdTcJ9PAXcHvJ+WE1Yvx4+NmH9u5P7zTIAirftQwkLwRwTiRIaRHCXC8n5XD9n3B60/EIAtFJHBqAE2UH6n5W7KzmH+pHBJNYJxmJSuKJlXVYeTA2o+08OZxXQLiGSIIsJs6LQqX40y/enABBqs5TZOd1qlXOFfXQcT6e+I5zE6HAsZhNXtG/Myt3HfGwW/3aIvhc2JCEyFICo0CBdR1ou7RrEEmQ1YzYpdGocz/Ld+o9vz+aJmE36rdauQSypOUV+NSj1m6JlJiOzUkH14Nm+FnPrys7dcvFg1F0bkAW5ZcdkfrburAFcxWhpJ1AiogEw97hCb6mXdhfIwtxKafroqNsELSsVmZMOmoq6eyPm5p0q2Zm7XI5nz2ZkYV6FBExgtoJQwGJFFmSX6vLWC2GxcRZma1dGU2v+Oc+otoUthGiJvspUPfTSTQIWSin3/BWCUtMziY91lH2Pi3WQmp5BWnoG8bExPsd37t4HQGZ2DjGOKABiHFFk5VRTGSMdaJlpZd+1zHTMTVtVsrN07oX9xlGIiDoUvvBQtXGViChkjv7YL3OyEOGRldIs0xAVg5aZ7k0nKx1TM/+rL1r7Xon7t0162lkZFH/3GRFvLkC6SnBv34xnxxY93br1MbdqS9BNdyLdLpzvv4l6aJ/fdM3x0XiSM8q+u1MyCGrfwtcmzoE72avVk5KBOc6BVFXUrFwSnp+AvWVjincdJOXp2UhnCQAxE/5DxNB+qPmFHBvufzngtIJi4sOCyr7HhdnZmVz5h2ZHUjbD3l1NTKidCZe1oqkjjKaOUF5bs48cpwub2cTaw2m0jve2Hl9ds5fvd58k1GZm7g3d/WoASMstIj4ixKsjIpidx9N9bI5m5OFRNe6c/QNFJW5u7tWGqzo1o2lcJK/9uIWcwmJsFjNr9x6ndWJMxSz4ZvN+rmjf2K8GERGNzMks+y5zM1EaNPO1CY/C3KYbzjlPYLuuacUkdJvIGJR6jVCP7QdAcdTF1KgV1ituBo+bkh/eQztx0L+OsEhknrcLS+ZnodRtUsnG3KIzxR9Ox1q3cTnbbNwbFhF83yvgdqH+sRP1sHc5actl12Nu1wuKi3B+ON2vhtPm/2sLWwgxFX1ZQAFsQp9uKYBPAi3ILYS4SwixRQixZd77n/gzqxJZRUtMCFFlA02c7rpeVUasnIF7y1ryJ99G4UuPYb/+jlOKe1r4aYWa23TA1ndwWV+0CAnF0uVicu+5kdy7rkXYgrBeMkAPU0yIkDDyHx6L84PZhEx8sppMqzifCjqqPmWJMJmwt2lK9seL+GPIvWjOYhyjh5WZpL/8PgcvvY28hSuJvPUqvwqqOuuKWbaKC2fx6L4suP1SbrzoAiZ8rf9ANY4OY0TXxty94Bfu+WITzWPDMSnean3vJS1Zcnc/Breqx6e/HiUQVevwVaJqGntOZvDaHZfzxsiBzFm6jaPpuTSOq8OIy9px99wfuWf+jzRPiMak+Madu2wbJkVhcEdfx1ctFYTZrrqDksUf+HdMVjv2Wx+gZOHbUFLaolVMEBSK8/UHKfnhPey3TAqcZw1uLuuAW3Et/7RyvbUHY25xEUWvTaDolXvBYsN0offp0b3yc5yzxuPZtR5L5wHV5nPK/L1T0/9Wqmth3wm0kVK6yx8UQswAdgPPVRWp/ApYp7q8anysg5Q0b4svNS2DWEc0bo+HlLR0n+MxDv1xLzqyDukZWcQ4okjPyCKqTuD+OS0rHSU6tuy7Eh2Dlp3p117duwMlri4iLDxgXC03C1FHb2WLOlHIvOxqNHhbYEpUDDIro5KdqUFjgu+eQsH0qcgC/bHT3LYTWloyMk9/knD/shpTizaw5me0rHTcv6zRdR/cC5qGCI8os62IJyUDc4L3icYS76j0ctCdkoElIYY/H2bN8Q48aZlIqYcVb9db8Hk/rsUx+vpKeeR+t5L6c58kY9ZHVWqIC7WTku99VE7NLyamwsvB8t0blzSOZfrPu8guchEZbGVouwYMbdcAgFmr9xIXVvnF4qBWdbn3q82M7dW8Sg2gt6hTcgu9OnKLiAkPrmATQp1gO0FWC0FWC50ax7MvOYuGMREM7dqCoV31p5NZi7cQF+GNu3DLAdbsOcZbdw1GBHCGMjcTUSe67LuIiPZp6QIoiU2w3zRRDw8Jw9SyEyWqivr7JlBM2IdPwbNtNeruX3zSVXfpXYjaiYO6kw0Jh4pdGX/a52UhwqO8OsKikPm+9Vmp2wjb0HF6eHAY5qbtKdE0MJnQctKhKB8Ade8WTInNUHet84nv2b0e+w2Tca/+ym95nBa1+KVjdX3YGlC3iuMJpWFnnct6dWfhj8uQUrJ91x5CQ0OIcURxYcvmHDuRxImkFNxuN4uXraJPr+5lcb5dvBSAbxcvpc8lPQJlgXpoL0p8PZSYeDCZsfboi3vreh8bJc572qYLmiHMZmR+XsC47q3rsV56BQDWS6+olKaPhoP7UBISUWLjwWzGcnFfXFt87YUjlpAp0yh8dTpasvelk5aRhrlZa7DqL/jMbS8qe1np2rQWc9uO+jkkJCLMFr/OGsC5cz/WC+piSYwDi5nwKy8lf5nv+swFy34h4pp+ANg7tEDLL8STno2akY0nOR1ro3oAhPToQMlBvc/X0tBbfmH9uuE6XOGlWTnaJERwLLuQkzlFuFWNJXuT6N00zscmo6C47OlrZ3IOUkrqBOlOPKtQ74JJznOy/EAKg1rpeo5me53vqkOpNIoK9asBoE1iDMcy8jiZlY/bo7Jk+2F6t27gY3NZ64b8diQFj6rhdHnYeSyNxrF6A+HPF5DJ2QUs33WEQR30lvS6fSd4d+UOZt4+gCBr4DaSduIgSnQCIjIWTGbM7Xuh7tnsY1P0/BiKnr+boufvxrNzAyXfzNGdNWC77h60tJO41/iO6PDs/gVTE/2Fq3AkgMns11kDaEmHUaLiEXViQDFhatMdz/5ffWycr03E+doEnK9NwLNnEyWL30XdvxWZm4mpXlO9DxtQGrVByzip5x3pva6mZhehZSYHLI/TQUq1xp/zjepa2PcDy4QQB/Auxt0AaAqMO50MpzzxHJt/20FOTh79rrmVsXcOx+PRH01uGHoll/bowpoNmxk07A6C7HamPayPRjCbTTw8YQyjJz6KqqoM/dflNG3cEICRw4cx6bHpfPX9EhLiYpjx9COBRWgazndnEfLQC6AouFYuRjtxBGt//bHdtfQ7LF0v1Z2vx4N0lVA466mAcQFKFn5C8PgnsF42GC0zjaKZTwbQoFI0/xVCH3lRT2dFqYYBV+safl5I0HW3IULDCR6llwGqSv6Do1EP7sG1cRXhL8wFVcVz5AAlS7/X461YRPCYqYS/9A7S46bw9WcDl4WqkfLfN6n/9tP6sL4vfsJ18Bh1bhoMQM4niyhYuZmQ3l1osmy+PqzvwZfLoqdMm03dlx5AWMy4j6eQVBoWO2WE7sg1iTspze8IEQCzovBg/wsZ88UmNE0ypG0iTR1hfL5N/xG6vkNDlu5PYcG2o5gVgc1s4rmrOpa1VCd9u5XcYjdmRfBQ/wsJt+uOfNaqvRzJLkBBkBARxCMD/I8QATCbFB4c0oMx837UdXRpTtP4SD7foL+uub5HKxrH1aFn80SGvfw1QsDQri1oGq+3RCe9v4zcohLMJoWHrulJeLD+g/rcN+txeTTunvsjoL94fPRaPy+kNY2Sb+cRdOfjoCi4Ny9DSz2OudvlAHh++cmvfuWCllg6XYaafISg8S8B4PrxI9R9v+LZshzbdfcQNGEmqB5KFswKWBZIDdeP72G/6QFQFDzbViEzTmK+qK+u49flfqNqSYfw7NlE0MinQVPRUo/i+W0FANa+N6BEJ+jD/XIzzv4IkVLttRVRVZ+xj4EQCvp2N/XQuxZPAJtlDX+ejB1nvBg7zngxdpzxYuw44+Vs7DhT/OvCGt9o9ouuPuP8/k6qHSUi9bUKz32NMjAwMKgJtbiFbUxNNzAwqF2o7uptzlMMh21gYFC7qMWjRAyHbWBgULswukQMDAwMzhNqcQvb2ITXwMCgdqFpNf9UgxBioBBinxDiYFWzu4UQU4QQ20o/u4QQqhAiqjTsiBBiZ2nYlrNxakYL28DAoFYhz9JLRyGECXgdGEDpcGYhxEIp5e9leUn5IvBiqf1VwAQpZfmpqX2klJWnMJ8mRgvbwMCgdiG1mn8C0xU4KKU8LKV0oa+rNCSA/U3AqS2edIr85S3sf8KklZBX559rCQAUTRh1riVQ/7qg6o3+BnLe/e1cS0B1/zPmTGSkBJ4y/3dgt/nfXOLvpOWjZyGRs9eHXQ/vDG/QW9ndqjIUQgQDA/GdAS6Bn4QQEnirdI2lM8LoEjEwMKhdnMIoESHEXcBd5Q7NKedYT2VpzquAdRW6Qy6WUiYJIWKBn4UQe6WUq2ssrgoMh21gYFC7OIUWdvmVRavgBFC/3PdE9P0AquJGKnSHSCmTSv+mCSG+Ru9iOSOHbfRhGxgY1C7OXh/2ZqCZEKKREMKK7pQXVjQSQkQAvYFvyx0LEUKE/fk/cDmwq2LcU8VoYRsYGNQuPGdnYwIppUcIMQ5YApiAt6WUu4UQd5eG/7mL+FDgJyllYbnoccDXpStKmoGPpZQ/nqkmw2EbGBjULs7iTEcp5SJgUYVjsyt8fxd4t8Kxw0D7syakFMNhGxgY1C5q8UxHw2EbGBjULoy1RAwMDAzOE4wWtoGBgcF5gtHCPosZtu9C0H/GgWLCteIHShb6zuQ0d7qYoGEjQJNITcX5/muo+3YFjCtCwgge/ziKIx4tI4WiV/6LLCzwq+HR6TNYvW4TUZF1+ObD2ZXCpZQ8O3M2azZsxm638cwjk2jdoikAazdu4bmZs1E1jWuvGsjI4cMAyM3LZ9Jjz5KUkkrd+DhemvYQEeFhgcuibRfsw+/R9+5buYiS7z/1Db+oJ/ZrR4DUkKpK8UdvoO7XyyJo5GTMHbsj83IoeGhkWZygex7FlKAPHRXBociiAgoeHR1Qh6l5B2z/uqNsD0H3qq+rtFMSmxA05lmKP5mh78BtthB01zQwW0Axoe7agGvpZwDYbpqI4tA34hVBIUhnIc5XJ/vVYOvWhfDx+rUt+v4HCj/0rRemBvWp8/BULM2bkT93PoWfLCgLCxl2HUFXXQlS4jl8mJzpz4PLjb1Pb0LvuB1zwwZkjhqDe9/+gOUAYOvehToTxyEUhcKFi8h/v0L9bFifyMcewNqiGbmz36bgI6+O+K8/RhYVITUNVJW028cAEPX0Y5gb6tdECQ1FKyggbfhd+COs90XUe2IkwmQi89OfSHvzy0o29Z4cRXifzmjOEo5Nnolz12EsCQ4avHw/lphIpCbJ/HgJGe/4bsYbc9c11HvkDnZ2uAU1Oz9gWYRc0onYR0bre31+voSsOZ9Xsol9dDShvbuU7vU5g5LfD2FtVI+6M73rJFnqJ5Dxygdkv1c26o2oO/5N7IMjOdDtRtRs/5sBnxZnaZTIP5G/12ELhaAR4ymcPgUtM52wZ2bj3roe7eTRMhPPrq3kb10HgNKgMSH3PUH+5NsCxrUNuRnPrl8pWfgJtqtvwnb1zRR/4n8W6DWDB3DztVfz8LT/VRm+ZsNmjp1IYtFn89mxey/T/vcan8ydiaqqPP3S68ydOZ34WAc3jBxPn17daNKoIfM+WED3zh0YOXwY8z5YwPwPFzBxbIBp+ULBftt9FD7/ADIrndCn3sD96wa0pHJlsftXCn7Vd1JX6jcmeNxjFEwdAYBrzRJKfv6W4Lun+iTrfP3psv/tN92NdBYSEKFgu3oUzvlPIfMyCbrneTx7NiPTTlSysw4cjnpgu/eYx41z3pPgKgbFRNDdT6Ps+xXt+AFKPplRZmYdfBuyuMi/BkUhfOJ4siZMQU1LxzFvNiVr1+M54i0LmZdP3sxXsV/ayzeqw0Hwdf8m/dbbweWizlNPENSvL87FS/Ac/oPshx8n4oGJgcugnI7IKeNJv1fXEfvumzjXrMfzh1eHlpdPzkuvEdS76k1008dORMv1dUBZj04r+z/ivrvRCgNcE0UhcdpoDt3yOO6UTJovfIncpZsoOeCdIR3WpxO2RnXZ03s0wR1bkPj0GA5cMwWpqiQ9/TbOXYdRQoJo/v0M8tduK4trSXAQ1qsDrhNpNSqLuCfGcnzEI7hTMrjgy5kULNuI65BXR0jvzlgvqMfhASOxt29B/H/HcfT6Cbj+OMmRIfeWpdN0zfvk/7yhLJ453kHwxR1xn6yBjtOhmn1qz2f+1okzpqYt0VKS0NKSQfXg2rAcS+cKFb+kuOxfYbPz50zQQHEtnXriWr0EANfqJZXTrEDnDm0Dtn5XrN3I1QP7IYSg/YWtyM8vID0ji5179tMgsS716yVgsVgY1K83y9fo212uWLOBIYP6AzBkUH+Wr97gN30AU5OWaKknken6+bg3rsDSqWfgsihXEdV9O5GFgVsmlm69cW/wv7s1gFK/KVpmCjI7FVQPnu1rMbfqUjmtnoNQd21EFuT6BrhKNZpMoFT9+29u2xPP9rX+dbZqiXoiCTUpGTwenEuXY+vlew21nBzce/chq2g9CZMJYbOBSUHYbKgZmQB4jh5DPX68kr0/rK1b4jlx0qvj5+UEXep7TbTsHNx79oGnRntQVyKo/2U4f/J/TYI7NKPkSDKu46lIt4fs79YQMcB3+YqIAd3I+lLfhbzot32YwkMwx0biScvGueuwrrPQScnBE1jiosvi1Xv8TpKefbdGDs3erjmuo0m4j6eA20PeD6sJ7d/Dxya0X3dyv14GQPH2fShhIZhiIn3Pp0d7XMdS8CR5nXPsw3eR/uLbf51jPYvLq/7T+Ftb2EqkAy3Te+G0zHTMTVtVsrN07oX9xlGIiDoUvvBQtXGViChkjj6FX+ZkIcIjK6V5KqSmZxIf6yj7HhfrIDU9g7T0DOJjY3yO79y9D4DM7BxiHFEAxDiiyMqp4NgqICIdyKx07/lkpWNqUrkszJ0uxj5sJCK8DkUvPVLjczC1aIuWm42WejKwjvAoZK539UeZl4VSv1klG3PrbjjnPYktsWmFBBSCxr2AEh2Pe+OPaMcP+AQrF7RGFuQgM5P9a41xoKaVu7bp6VhaVy6LqtAyMij4dAGxX36GLCnBtXkLrs2nt/SwKdaBmurVoaZlYG1TMx06EsesFwFJ4dffUfjNDz6h1g7t0LKy8Rz3f00s8dG4k73Xw52cQXDHFpVtkrx1x52SiSUuGk9atjevxFiC2jSmaJteP8P7d8WdkknxniM1OhNLXDSeFK8OT0oGQe0r6Ihz4Enx6vCkZmCJc6Cme3WEX9mbvB9Wln0P7dsNT2omJXv/qJGO0+I8dMQ15bRb2EKIEQHC7hJCbBFCbHn3YFL5gCqsK//KuresJX/ybRS+9Bj26+84pbhnA1nFL78QosoGQZWyakKVp1M5A8/WdRRMHUHRzMexX3t7jZO39OiLe+OK0xNSQYftXyMo+fGDql/mSA3nq5MpfO4ulMRmKHH1fYIt7XsFbF3rEqrX4DdqWCj2Xj1JH3YTaddch7DbCbq8f43iVpHaaesASBt1H2m3jSbj/gcJue4arB3a+YQHX96XogCt6xprqKbuKMF2Lpj9ICefmodW4ETYrcSNu57kGR/X7ESgZtekSpNyNhYzof26kb9Yv/7CbiN6zI1kvPJBzXWcDmdvavo/jjPpEvmvvwAp5RwpZWcpZefbm9YtO65lpaNEx3ozj45By870m4G6dwdKXF1EWHjAuFpuFqKO3roVdaKQedlVpldT4mMdpKR5WxepaRnEOqKJi3WQkpbuczzGoT9yRkfWIT1Db+WnZ2QRVSciYB4yKwMR5W2tK1ExyJwAZbFvp14WoeHVn4CiYOl8SY0ctszLRER4nyZEeBQyz3epTaVeE+w3TST4gTcxX9gd25C7MLXu6ptQcRHqH7swNe/oo8PUphueHesCalDT0jHFlru2MTFl3RrVYevcCTU5BS0nF1SV4tVrsLS9sEZxq9QR59VhinWgZtR87XmtVLOWnUPxyrVY27T0BpoUgvr0wrk08DVxp2RgSfBeD0uCA3eq7/VwJ2diqeutO5b4aNxppTZmExfMfpDsb1aR+6PeLWdrmIC1fhwtF79C67VzsSQ4aPHDTMwxdQLqMMd7dZjjHd48fGy8OsxxDjxp3usWemlnSnYfQs3MAcDaIAFLYhyNFr5Ok+XvYI53cMHXszA5zuyJuBKqWvPPeUZAhy2E2OHnsxN9rvwpoR7ai1o3l/oAACAASURBVBJfDyUmHkxmrD364t663ldQnNfBmy5ohjCbkfl5AeO6t67HeukVAFgvvaJSmqfKZb26s/DHZUgp2b5rD6GhIcQ4oriwZXOOnUjiRFIKbrebxctW0adX97I43y5eCsC3i5fS55IegbJAPbwXU3w9ROn5WLr3wf1rhbKI9ZaF0rAZmCzIgurfqJvbdEJLPobMrt7ZaCcOojgSEJGxYDJjbt8LdY9vl0LRi2MpemEMRS+MwbNrIyXfzkH9fROEhIM9uDRTK+Ym7dDSvY/7pqbtkOknK/0AVMS9dy+m+vUwJcSD2UxQ/76UrKvZNVRT07C0aQ02GwDWThf5vKw8FVx79mIur2NAX5zVvIv4E2G3I4KDyv63deuM+5D3sd/WpROeI8dR0wJfk6LtB7A1qou1fhzCYibyqkvI+/kXH5u8pZuIurYPAMEdW6DmF5V1hzR44V5KDp4gfZ53REbxvqPs7vQffu81it97jcKdnMG+K+/Hk57jV0fxzv1YL6iLJTEOLGbCr7yUgmUbfWwKlv9CxNB+ANjbt0ArKPTtDvlXb/K+X1X2vWT/EQ72uJlDfUdwqO8IPCkZHBl6H2rGmTWwKvH/uA87DrgCqFiiAjh1r6hpON+dRchDL4Ci4Fq5GO3EEaz9rwLAtfQ7LF0v1Z2vx4N0lVA466mAcQFKFn5C8PgnsF42GC0zjaKZTwaUMeWJ59j82w5ycvLod82tjL1zOJ7Sl1k3DL2SS3t0Yc2GzQwadgdBdjvTHp4AgNls4uEJYxg98VFUVWXovy6naeOGAIwcPoxJj03nq++XkBAXw4ynq+lv1jSc779KyJTn9eF0qxejnTyKte+/9LJY/j3mLpdi7TUAVA/S5aLode9og6Cxj2Bu1R4RGkHYK59S/NV7uFctBsDSo0+1LxvL6yhZOI+gOx4DoeDeshwt7TjmrpcD4Nn0k9+oSlgktuvHgTCBEHh2rkfdu7Us3NyuF+7qukMAVI28GbOImqFfW+cPi/H8cYTgIXq9KPr2O5SoSBzz3kKEBIMmCbn+OtJvvR3373soXrGKmLfnIFUV9/4DFC38HgDbpb2IuP8+lDoRRL74LJ4Dh8ia9EBAHTn/exXHrOcRionC73QdIUN1HYVf6zpi35uNUqoj9MZrSb1xBEpEBNEv6HVVmEwULVlGycbNZUkHD+hTg+4QXcOJx9+i8ftPIkwKWQuWUnzgONG3DAQg86MfyVu+hbA+nWi1+q3SYX2zAAjp3Iqoa/vi3HOEFotmApD04gfkr9jqN7tAOlKfepP6858Gk0LuFz/hOniMOjcOBiDn00UUrtxMaO8uNF46H81ZQspDL5dFF3YbIT07kvLYq6ee95lyHjrimiKq6q8tCxRiPvCOlLLSXSeE+FhKeXN1GeTc1Oecj7ExdpzxYkr0/xj8d5K/5tzvcGLsOOPFbvtnjF1uuX/RGV8U57yJNfY5QSNn/DMqQQ0J2MKWUvodSFwTZ21gYGDwdyO1c95G/MswpqYbGBjULmpxl4jhsA0MDGoX5+Hoj5piOGwDA4PahdHCNjAwMDhPMBy2gYGBwXmCsfiTgYGBwXnCWZw4I4QYKITYJ4Q4KIR4sIrwy4QQuUKIbaWfx2sa93QwWtgGBga1i7M0rE8IYQJeBwYAJ4DNQoiFUsrfK5iukVL+6zTjnhJ/ucOWnnP/ePJPmLACEPzy3HMtgYIxd5xrCQCUFFjOtQSswf+MySJNBrnOtQQ86cXVG50vnL1RIl2Bg6U7oCOE+BQYAtTE6Z5JXL8YXSIGBga1CqlpNf6UX1m09FN+K6B6QPkF1U+UHqtIDyHEdiHEYiFEm1OMe0oYXSIGBga1i1PoEpFSzgH8bU9VkzWdfwUaSikLhBCDgW+AZjWMe8oYLWwDA4PaxdlbD/sEUH6B90QgqbyBlDJPSllQ+v8iwCKEcNQk7ulgOGwDA4PahSZr/gnMZqCZEKKREMIK3AgsLG8ghIgXQt/tQQjRFd2nZtYk7ulgdIkYGBjULk5zv82KSCk9QohxwBLABLwtpdwthLi7NHw2cB0wRgjhAZzAjVJfArXKuGeqyXDYBgYGtYuzuPVXaTfHogrHZpf7/zXgtZrGPVMMh21gYFC7MJZXNTAwMDg/kMZaIgYGBgbnCUYL+yxm2KErwSPGgWKiZNkPlHzzsU+4tVd/bNfcpH8pdlI092XUo4cAsF15HbZ+V4IE9dhhCt94Htz6LDHbwKHYBg0FVcX960acH74VWEfbLtiH36Pvp7hyESXff+obflFP7NeOAKkhVZXij95A3b8LgKCRkzF37I7My6HgoZFlcYLueRRTgj6SRwSHIosKKHh0tF8Nj06fwep1m4iKrMM3H86uFC6l5NmZs1mzYTN2u41nHplE6xZNAVi7cQvPzZyNqmlce9VARg4fBkBuXj6THnuWpJRU6sbH8dK0h4gIDwtcFn/BNQmZ8DhK3QY+ZZE/ZST+COrZmagHxoKiUPD1YnLf+cwn3HJBfaL/Oxlbq6Zkv/YOee9/AYApLgbH0w9gio4CqZH/5SLyP/5aj9O8MdGPjEcJDsKTlEL6w88hC4sCloWtWxci7h8HJhNF3/1AwQef+JZVw/rUeWQqlubNyHtrPoWfLNB1NKhP1FNly0hgqpdA/tx3KFzwJWGjRmC/5GLQJGpONjlPP1+2w3pVmFp3wj5sjF431/2Ia8kCXw3tu2O96ja9r1ZTKVnwFuqh0vdZQSHYh9+PUvcCkJLi919G+2OPXh6XXY31squRmoq6axMlXwXeOs/coSvBd5SrF19XqBeX9Mc2tLReOJ0UzdHrhVK3PiETn/CeT1wCzk/foeQH/ZrZBpXeq5qKe+tGnB8EvldPGcNhnyUUheA7x1MwbTJaVjphz87GvWUd2gnvLtdqWjIFT4xHFhboFWb0JPIfHouIcmAbfC15E24Dl4uQCU9gvbgvrpU/Ym7TAUuXXuRNuhM8bkR4NfsWCgX7bfdR+PwDyKx0Qp96A/evG9CSvDo8u3+loHQXc6V+Y4LHPUbB1BEAuNYsoeTnbwm+e6pPss7Xny77337T3UhnYUAZ1wwewM3XXs3D0/5XZfiaDZs5diKJRZ/NZ8fuvUz732t8Mncmqqry9EuvM3fmdOJjHdwwcjx9enWjSaOGzPtgAd07d2Dk8GHM+2AB8z9cwMSxfnd6+8uuSeHLT5XFD/rPGGRRgLJQFKIeupfUu6fiSc2g7kevUbRqA+7Dx7wacvPJeuF1gvtc7BtXVcl+6S1cew8igoOo+8kbFG/civvwMRxPTCRrxhxKtu4gdMgVRNx2PTlvvBdQR8Tk8WSOn4Kalk7M/NkUr1nvswu7lpdP7suvYr+0l6+MY8dJv31UWTpx335O8Wp9K9SCjz4jf+47AIRc/2/CRvyH3BdfpkqEgv2meyh65WFkdgbBD83Cs2MjWrK3LDx7t+HZru9grtRrhH3UwxQ9qedtH3Y36u6tFM95BkxmsOq7yZuat8PcvgeFT4/R75GwCP/lUHoOwaPGU/DUZLTMdMKen417cxX14rHSetGxK8F3TyL/obFoScfJnzzSW6ZzvsC9aQ0A5gs7YOnai7yJNbxXT4davIHB3zoO29S0JVrKSbS0ZPB4cK9bjrWz7w2o7t+NLCzQ/z/wO0p0TFmYUEwIqw0UE9jsaFkZANguH0LxNx+Dxw2AzMsJrKNJS7TUk8j0ZFA9uDeuwNKpp69RiXdtBWGz+yzZqO7biSzMC5iHpVvvancu79yhbcDW74q1G7l6YD+EELS/sBX5+QWkZ2Sxc89+GiTWpX69BCwWC4P69Wb5Gv0GXrFmA0MG9QdgyKD+LF+9IaCGv+qalMfaow+utcv8arBd2ALP8SQ8J1PA46FwyUqCL/O9Hlp2Dq7d+8Hju/6HmpGFa+9BAGSRE/fhY5hiHQBYGiZSsnUHAM6NvxLc75KAZWFp3RLPiSTUJL0snEuX6y3jCjrce/ZV0uFzPp0vQj2ZhJqSWqrL26oXdt+6VBHlghZoacnIjBRQPXg2r8LcroevUbm6ibVcevZgTM3a4l73o/5d9UBpo8HS+196S/3PeyQ/N2BZlNWL1NJ6sXY51i4V6sW+cvViv2+9+BNz24vQUk+ipetlYbtiCMVf1/xePR2kJmv8Od+otoUthGiJPgf+lz9n9JQeHyil/PFUMlOiYtAy08u+a1npmJq19mtv7Xsl7t82ASCzMij+7jMi3lyAdJXg3r4Zz44terp162Nu1Zagm+5Eul04338T9dA+/+cU6UBmVdDRpFUlO3Oni7EPG4kIr0PRS4/U+DxNLdqi5WajpZ6scZyqSE3PJL7U+QDExTpITc8gLT2D+NgYn+M7d+vnm5mdQ4wjCoAYRxRZOYFvzL/qmvyJuVU7vSxS/JeFKdaBJ8WrwZOaga1ty4C6q8JcNw5ry6aU7NwLgOvQEYIu64Fz5QZCBlyKOb6yQ/HREeNATU0r+66mp2NtXbleVEdQ/74U/ez7AxU2+k6CB16OVlhI5rgJfuMqkdFo2eWuR04GpkYtKtmZO/TEes0IlLA6FL2md8UojnhkQS722yah1GuEeuwgJQveBFcJSmw9TE3bYBtyG9LtouTLeWhH9/vXERWDlnEK9aKft174HL+4L6613oaLklDhXn0v8L16WpyHjrimBGxhCyHuA74F7gV2CSGGlAueHiBe2YIq7x6uZjamn9aGuU0HbH0Hl/VFi5BQLF0uJveeG8m961qELQjrJQP0MMWECAkj/+GxOD+YTcjEJwPnWeUs/8o6PFvXUTB1BEUzH8d+7e2B0yyHpUdf3BtX1NjeH7IKTUKIKotMVHVOp59xlYdP5Zr8ibVXv4Ctaz2hKsSf4iL0IshOzP8eJ+vFN8v6qTOfeInwG4aQ8PHriJAgpLu61fkq66jqGgTEbMbWqyfFy1f5HM5/az6pQ2/AuWQpIdcOPSUNVdbNbespenIUzjf/i+3q/+gHFRNK/aa4Vn1P0fRx4CrGesUNZWEiOIyi5++n5Kt5BI16OPB51PAeAb2bw9ZvcOW+aLMZS5eLca1f6U3WVHqvPjQW5/uzCZn0ZGAdp8NZXA/7n0Z1XSKjgE5SymuAy4DHhBDjS8P8uggp5RwpZWcpZefbG9ctO65lpfs8NilRMcgqHqFNDRoTfPcUCl54BFmgdz2Y23bSHxXzcvUXi7+sxtSiTVm67l/0PjL14F7QNES4/z46mZWBiKqgI8f/SyB1306UuLqI0HC/Nt7EFCydLzkrDjs+1kFKmrd8UtMyiHVEExfrICUt3ed4jCMagOjIOqRnZAGQnpFFVJ3AfZV/1TXREzNh6XoJrvWBy0JNTfdp/ZrjHKjp/q9HJcwmYl96gsJFyylavrbssPvIcVLHPEjyzfdQuHgFnhOBGw9qejqmuFjvOcfEBHw5WBX2Ht1w79+Plp1dZbjz52XY+1zqN76WnYESWe561HEgc7L8az64CyUmARESjszJQOZkoB3RW6yeX9dgaqC/pJY5GXi2rdPzOLIfpIYI9V83tMx0FEcN6kXDxgSPmULBc9568SeWjt1QD+9H5nrLQsuscK/KwPfqaXH2pqb/46jOYZvKLWxyBN1pDxJCzCCAw/aHenAfSkIiSmy8/ut7cV9cW9b72AhHLCFTplH46nS05BNlx7WMNMzNWpe9RDG3vajsBYhr01rMbTvqJ5SQiDBbdCfiT8fhvZji6yFi4sFkxtK9D+5ffXUosd4fGqVhMzBZKlXIqjC36YSWfAyZXblynyqX9erOwh+XIaVk+649hIaGEOOI4sKWzTl2IokTSSm43W4WL1tFn17dy+J8u3gpAN8uXkqfS3oEyuIvuyYA5nadUJOO+XQ/VUXJ7n2YG9TDXFfXEHLFZRStCtz3Xh7HE5Nw/3GMvA+/9DmuRJa+0BKCOqNuIf/z7wOm496zF3NiPUwJuo6g/n0pXrs+YJyKBA3oi/Nn33cXpkTvqpr2Xj3xHD1WMVoZ2tF9KLF1EdFxYDJj7tIbz46NPjYiJqHsf6V+UzCbkYV5yLxstKx0RFyinm/LjmUvKz3b1mNq0V6PH1uvtD4HuEcq1oteAerFLN968SdVPV2d6r16WtRihy0CPfIJIZYDE6WU28odMwNvA7dIKU3VZZB9/WU+GZg7diP49nGgKLhWLKb4qw+xDrgaANfPCwm+ewqWbpeiZegvKVBV8h/Uh8bZh92OtWdfUFU8Rw5Q9OaL+ssLs5ngMVMxX9AU6XHj/OBNPLt+K8tTsVeWaW7fFfstpcP6Vi+mZOHHWPvqm0a4ln+P9cobsfYaAKoH6XJR/Olb3mF9Yx/B3Ko9IjQCmZdN8Vfv4V61WA+76wHUg7/jWl7ZOVTcwGDKE8+x+bcd5OTkER1Vh7F3DsdT+jLrhqFXIqXkmRlvsHbjFoLsdqY9PIELWzUHYPX6TTw/aw6qqjL0X5cz+jZ9eFVObh6THptOcmo6CXExzHj6EZ8Xm1VtYPCXXBMg+J4H8ez/HdfPlde8yT3gu4FBUK+uRE3Rh7IVfLuE3HkfE3adfj3yv/geU3QkCR+/jhISDFKiFTk5+e+RWJs1IuHdmbj2Hy57ZM9+9W2cazcRdvNQwm/Qz6No2VqyZ/kOY6tqAwNbj25EjL8HTApF3y+m4L2PCL7mKj2Nb75DiYok5u23ECHBoEmk00nazbcji4oQNhtx33xG6nW3IAu9o2Iin/kv5ob1QdNQU1LJeeFltAzvD3pYO5uPBtOFXbBfP1qvm+t/wrX4UyyXDAbAvWYR1suvx9y9v/5SsbQ/+s9hfUpiY+zD7weTBS0jmeL3Z0BRAZjM2P8zESWxMageSr6ci7pve1meVW1gYL6oW+lwTwXX8sUUf/kh1stL68VPCwkeMwVL90vLXiiiquRPLR3GarURMWcBuWNvhvIjhMxmgsdOxdyo9F59z/dejfxy5Rl37uWNurzGnjh87k9nszPxL6c6h50IeKSUKVWEXSylXFddBhUd9rmgKod9LjB2nPFS0WGfC/4pO85UdNjngn/KjjNnxWHfOaDmDnv+z+eVww44SkRKWfk5xxtWrbM2MDAw+Ls5H4fr1RRjarqBgUHtwnDYBgYGBucJ599ovRpjOGwDA4NahfTUXo9tOGwDA4PaRe3114bDNjAwqF3U5peOxia8BgYGtQvtFD7VIIQYKITYJ4Q4KIR4sIrwW4QQO0o/64UQ7cuFHRFC7BRCbBNCbKkY93QwWtgGBga1irPVwhZCmIDXgQHACWCzEGKhlPL3cmZ/AL2llNlCiEHAHKBbufA+Usozn/Zcyl/usFO2Bf/VWVRL/euCzrUE4J8xaSX0zbfPtQQAVl346LmWQPBZ3Kz1TIhPCrxu+t9BQUnIuZYAQM/qTarn7F3WrsBBKeVhACHEp8AQoMxhSynLz9ffCCSetdyrwOgSMTAwqFVIT80/5VcWLf3cVS6pesDxct9PlB7zx53A4vJSgJ+EEFsrpHvaGF0iBgYGtYpTeXCSUs5B78aoiioXma3SUIg+6A67/FZEF0spk4QQscDPQoi9UsrVNVdXGaOFbWBgULs4ey8dTwD1y31PBCqt0SuEaAfMA4ZIKcvW45VSJpX+TQO+Ru9iOSMMh21gYFCrkFrNP9WwGWgmhGgkhLACNwI+S08KIRoAXwHDpZT7yx0PEUKE/fk/cDmw60zPzegSMTAwqFWcrXfJUkqPEGIcsAQwAW9LKXcLIe4uDZ8NPA5EA28Ifeckj5SyMxAHfF16zAx8fKpbKlaF4bANDAxqFVI9eyumSikXAYsqHJtd7v+RwMgq4h0G2lc8fqYYDtvAwKBW8Q8ZrfmXYDhsAwODWoXUzqs9CU4Jw2EbGBjUKowW9lkk5JJOxD06GmFSyFmwhMw5n1eyiXtsNKG9u6A5S0ieOoPi3w8BoISFkDB9PLZmDQFJ8oMzcW7bS8z9wwnt1x2khiczl+SpM/Ck+d9pGsDUvAO2f92h75u3eRnuVV9XaackNiFozLMUfzIDdddGMFsIumsamC2gmFB3bcC19DMAbDdNRHHom/eKoBCksxDnq5P9ajB36Fq6Z56JkmU/UPLNxz7h1l79sV2j79VIsZOiuS+jHtXLwnblddj6XQkS1GOHKXzjeXC7CJnwOErdBrqG4FBkUQH5Uyp1sZXx6PQZrF63iajIOnzz4exK4VJKnp05mzUbNmO323jmkUm0bqHvxL124xaemzkbVdO49qqBjBw+DIDcvHwmPfYsSSmp1I2P46VpD/nsK+mPtk//h7h+HVCdLn4dP5vcnUcq2TS643KajBpIaKN4FrUejSsrH4D4KzrRaur1oGloqsbOxz4ga5O+e3jjkQO54NY+IARHP1zOobmB3/20eOY2Yvp1RHWWsOu+N8mvQkdQgxjavTUec50Q8nceYec9ryHdKuawINq+MQ57PQfCpHDkze9J+nQVAA1GDSLx1r4AnPhoOcfmLK6ULkBo74uo9/goMClkffYz6W9+Ucmm7hN3EdanE5qzhBOTX8G5+xDCZqHJZ88hbBaEyUTu4nWkvqzXqbiJtxA+oBtIiScjl+OTZ1Z7jwA0mnYHdfpdhOZ0cfD+Vync+UclG1v9WJrPnoC5ThiFOw9z4N5ZSLeHyCu60OCBm0DTkKrKH4+/Q/6mvQDU6dOBRk/dASaFtI+XcfK1qu+/00HK2tvC/nuH9SkK8U+O5fjIxzk06G7C/9Uba9P6PiYhvTtjbViPQ/1HkvzYLOKfGlcWFvfoaApXb+XwwNEcvmocJYf0SUiZ877gj6vu4Y+r76VgxSYc424OrEMo2K4ehfOdZyh6+X7M7XshYquYUSoUrAOHox7wblaKx41z3pM4Z03COWsSpuYdUOo3A6Dkkxk4X52M89XJeHZtxLP7l4BlEXzneAqemUrehNuwXtwXJbGhj4malkzBE+PJn3wnzi/eJ3j0JF1WlAPb4GvJe3A0eZNGgKJgvVh3BIUvP0X+lJHkTxmJ+5dVuH8JPE7/msEDmD3jab/hazZs5tiJJBZ9Np8nH7iPaf97Tdemqjz90uu8+dI0Fn70FouWruTQH/qO6fM+WED3zh1Y9Nl8unfuwPwPFwTUABDXrwOhjeNZ2mMi2ybPo/3zVU/jz9q0j/XDplN03Hcn9vQ1u1jR90FW9H+Y3+5/i44vjQIgrGUiF9zah1WDHmNF3weJG3ARIY3i/epw9OtASKME1na/n98nz6X1C1X/2DV79GaOvvUD63pMwJ1TQL2b9fKvf8cVFOw7yYa+U9n876do8eRwhMVEaMtEEm/ty8aBj7Ch71RiBlxEcFU6FIV6T93NH7c/yf4B91Dn6kuxVbhHwi7rhLVRXfZdNpqTD79OvWfGACBL3By++REODLqP/YPvI6z3RQR3bKGXz5yvODDoPg4MHk/e8s3Ejb/Rbxn8SZ2+F2FvnMBvPcdxaMqbNH6u6sl6DR8dTtKc7/nt4nF4cguIvakfALlrdrK930S2D5jMwQlv0OSlsWXn2Hj6KH6/5Rm29b4fxzW9CGp+9mZ0n8Vhff84/laHHdSuOa6jSbiPp4DbQ94Pqwnr18PHJqx/d3K/WQZA8bZ9KGEhmGMiUUKDCO5yITmfL9EN3R60fH0NBq3AWRZfCbKX7Z7tD6V+U7TMFGR2KqgePNvXYm7VpZKdpecg1F0bkQW5vgGu0g1LTSZQqn5IMbftiWf7Wr8aTE1boqWcREtLBo8H97rlWDtf7GOj7t+NLCzQ/z/wO0p0TFmYUEwIqw0UE9jsaFmV15ex9uiDa+0yvxoAOndoG7D1u2LtRq4e2A8hBO0vbEV+fgHpGVns3LOfBol1qV8vAYvFwqB+vVm+ZqMeZ80GhgzqD8CQQf1ZvnpDQA2gt5CPLVgDQPavB7GEB2OLrVPJLnfXUYqOVz5Xtaik7H9TsLcOhDWrR9bWg6hOF1LVyNywh4TBnf3qiBnYmaTP9R+53K0HMYcHY61CR1SvNqR+p/8gJy1YTeyg0jSlxBxqB8AcYsedU4D0aIQ0q0fO1gNopTqy1+8hdnDlOhfcoRmuo8m4jqci3R5yvltN+OXdfGzCL+9OzlfLASj6bR+m0nsEQCvS66YwmxFmM39usu1zjwTbqr1HAKIGdiH9c/3poODXA5jDQ7BUURYRvS4k83v9GqctWEnUoK4+WirmGdqxKc4jKZQc088x49u1RF1RuSxOF00VNf6cb/ytXSLm+Gg8yd6bzZ2SQVD7Fr42cQ7cyd7WkyclA3OcA6mqqFm5JDw/AXvLxhTvOkjK07ORTv1GjZnwHyKG9kPNL+TY8EqrIPogwqOQuV4dMi+rrJVc3sbcuhvOeU9iS2xaIQGFoHEvoETH4974I9rxAz7BygWtkQU5yMxkvxqUqBi0TO95alnpmJq19mtv7Xsl7t826XqzMij+7jMi3lyAdJXg3r4Zzw7f1RvNrdqh5WajpZz0m2ZNSE3PJD7WUfY9LtZBanoGaekZxMfG+BzfuVvvgsjMziHGEQVAjCOKrJwKP3hVEJQQiTPJ+4henJxFUEIkJWk5NdaaMKgzrR++EZsjnA23vghA3t7jtH5wGJbIULRiF3H9OpCz/bDfNOwJURSfLJusRnFyFvaEKFzldFiiwvDkFSFVvYlWnKTbABybv4SOH0yh9443MYUGseOuV0BKCvYep+lDN2KJDEUtduHo34G8KnRY4qJxJ5W7R5IzCe7QvJKNq5yNKyUTS3w0nvRsUBSaff8y1oYJZH7wA85tZXM5iJs8nMh/90HLL+LQTQ9XW57W/2PvvOOjKNoH/p1rufTe6YSOgAYQkI5UG1jAhgJWrIgigu1VfAVRsYCKCHZR8bUA0qT3Lt3Qa0jvIbkkd7vz+2PDXY4klwQVJb/95nOf7O08z8yzs3vPzj4zOxMVQnGZcoqTM7FEh2IvUxemEH8cuQVQWhclyZl4RYU400MGdqTexLsxhwaQMPx1ALyiQig5XZngpwAAIABJREFUW8b+5Cz8rnT//f0ZanOnY5UtbCFERyFEh9LtlkKIsUKIQVXoOCdUmZd7umxKeeEL7vSiwrf3JcJoxNoqjuy5izlx0+OotiLCHhrqFEl/50uOdr+XvAVrCL77hqqOqko7vK4fSfHSryp+bpIqtunPUDDlQQx1mmCIdH9kNbft6rF1XSmVtHpMrdrh1XsQtq8/1qz39cPc4RpyH72d3AdvQXh5Y+nW103H0rVPla3r6plU3iYhRIWmVnjuqksFytVoBLqRvGQHK7s9w9aR07R4NnDuSBJHZizkmu8n0HnueHIPnEJ1KDXKt1wdVHj5aDJhvdqSv/8Ua9uMZnPv8bSYPBKjnzcFR5I4OWMB8fOeJ/7bCeQfOFXxUlYVVWI1ynfKqCpHBj1JQueR+LRtilfTek6R1Le+4mCXUWTPX0PYvdd7OOLzplTHlorOm0sma8k2dnd7gkOjpmrx7Ep0anyyPSBVUe3P5YbHFrYQ4mVgIGASQixHm+d1DfCcEOJKKeV/K9IrO6FKQpNBzjPhSMnAFO1qrZmjwsp1fNhTMjBHh3P+Ac4UFYYjLRMptbSiPVorLm/pBsIeuq1c2bkL11D3k/+Q8f43lR6XzMtEBLrsEAEhyDx3OwyxjbHeMVZL9/HH2OwqilUV5Y9tLqGiQpQT+zE2vRI1tXRSL4MBY6urKZkxrtLyQWtRlw1xGELCkRWENYz1GuHz8DjOvT4eeS5Pq5Mr4lHTkpF5WsvVvnUdxmatYP3y0syMmDt2wzb+IY82VIeoiDBS0lx2paZlEBEWit3hICUt3W1/eFgoAKHBQaRnZBEeFkJ6RhYhQYEV5t1wZF8a3NULgOzdx/GOcbXMrNEhFKVkX5TNmVsO4tsgAkuIPyVZ+Zz6dg2nvl0DQIsJwyhKznSTrzuyH7GlnYF5u49hjQ11s6P4AjvsmfmYAnwQRgNSUbHGuGRibu/Biena28u2k6nYTqfh2ySGvF3HODt3NWfnrgYgbuLtFCe52wGl139Mmd9IdCj2cr+RTCwxYRSWfrdEhWJPdZdR8wo4t2Uf/j3iKT582i0tZ/5aGn76srNDsixRIwYQeZcWzjq35yheMWHkl6Z5RYdSkuJejiMzD1OgLxgNoKhYokMpSS1/3vK2/IG1QSSmEH+tpR7rOkZLdAglqVV3gFaXv9D3/+uoqoV9K3AN0B14FBgspXwV6A8Mq2lhtn2HsTSIwVwnEswmAq7rTv7KLW4y51ZuJXCw1mlhbdcMNb8AR3o2SkY2juR0LA212Q19O7ej+Kh2IZrrxzj1/ftcTcnxRI92qIlHMYRFI4IjwGjC1LYrSoJ7SKHwzUconDqawqmjcezfQvH8WZqz9g0Aa+kc3yYLpsZtUNNdYQdjXBtk+tlyN4ALUY4ewhBdB0NEFJhMmK/pTcmOTW4yIiwC33GTKJj+Omqy65jUjDRMTVqCxUsz44qrUBNPOdNNbeJRkk4js9w75i6Gnl07sWDpSqSU7NmfgJ+fL+FhIbRu3pTTiUkkJqVgt9tZsnItvbp2curMX7ICgPlLVtCrW+cK8z7x2XJWXzuR1ddOJHnpDuoN7QZA8FVxOPJtNQqH+DaIdG4HXtEAg9nkHEFiCQsAwDs2lJhBHUj82T2mfuaz39jS5zm29HmOtCU7iLmtu5ZPfByO/EK3cMh5sjb+QeQNWmw5Zmh30pdq10/R2UxCu7XWyg0PxKdxDLZTaW52WGNDiRzUgeSfN5XLt3DPEedvRJhNBN3Qnbzl29xk8pZvJehm7Qbjc2UzlPxCHOnZGEMCMARo81oLLwv+17Sj+Jh23VgaRDv1A669mqJjFf9GUj5fyp6+z7Cn7zNkLdlG+G09APC7qgmO/EK3cMh5cjfuJ/R67RxHDO1J9lLNXmsDV6eq7xUNEWYTjqx8zu0+infDaLzqRiDMJsJu6krWsr9kQRbg/3ELG+29eAUoFEIck1LmAUgpbUKImvexKiopr3xE3U9f04b1/e83So6eJugOLcKS8+1izq3Zjm+PDjReOUcb1vfcO071lEkziXn7WYTZhP1MCkmlaRHjRmqOXJXYk9JIeWmGZztUleIFs/Ee9SIIA/Ydq1DTzmDq2E876G2/Vapq8A/G67bHQBhBCBz7NqEc3OlMN7Xpir064RBVoXDOe/g9/yYYDJSsXoKaeBJL3xsBKFm+AO9b70X4BeDzwFOl9aeQ/9xDKEcTKNmyloCpn4Ci4Dh5hOIVvzqztlzTm5INq6q2ARj38hS279pLTk4efQbfzSP3DcfhcAAwbMh1dO/cgfWbtzNw6Ci8rVYmTdRsMZmMTHxqNA+NfQFFURhyfT/iGmmjXO4fPpSnX3ydn35dRnRkONNee75KO1JX7CayTzv6bnkHh62YXWM+dqZ1+uZZdo+dRVFqDo3u60+TR6/HKyKIXqumkLpyN7uf/oSY6ztS97ZuSLsDpcjO9oemO/U7zh6DJcQPaVfYM+Ez7LmVLxiQsWIXYX3a0XXreyi2Yg486RrqeOU34/lj7CyKU7M58tpc2nz8BHHPDSNv30kSS1vOx6f9RKv3R9N5zVSEEByZNBd76Y2j7ZyxmIP9kA6FhAmfabHfC1FUkl6aSaMvXwGjgex5Kyg+cpqQuwYAkPXNUvJX78C/V3uarZ2lDesb9x4A5ogQ6r49BgwGhMFAzqIN5K/aDkD0+BF4NYpFqir2s+kkPv9Blecke+XvBPW5iqs2f4BiK+boUy6dFl8/z9GnP8Sems2p176m6cynqDf+Dgr2nyD1Wy0UF3pdJ8Jv64m0O1CLSjj88DTnMR6fOJuW376IMBpI/W4VtsNnKjLhoqjNw/pERTFKZ6IQW9GWuCkUQhik1AK6QohAYLWU8qqqCigbEvmn+LesOGM/Wv4R+FLzb1lxZpG+4oyTKO9/w4ozln/aBAC6JP/4p73t4RYDqu1zmiYsvay8e1Ut7O5SymKA8866FDNw799mlY6Ojs5FUptb2B4d9nlnXcH+DOAvW1hSR0dH56/icoxNVxd9LhEdHZ1aRW0eJaI7bB0dnVqF3sLW0dHRuUxQ1Nq78qHusHV0dGoVtTkkUntvRTo6Ov8vUaWo9qcqhBADhBCHhBBHhRDlJikSGu+Xpu8VQlxVXd2LQXfYOjo6tQopRbU/nhBCGIEP0KbnaAncIYS4cIa2gUCT0s+DwEc10K0xusPW0dGpVUhZ/U8VdASOSimPSylLgO+Amy6QuQn4UmpsAYKEENHV1K0xf3sMu8H45n93EVWS8/muf9oEAIrPmf9pE1j7L3jDEOC6/ZUvmnCpcKwpP/nRP0Hxz9WbRuDvxNyi4gm6LkeqE+o4jxDiQbSW8XlmlU5eBxALlH1nPhFtAryyVCQTW03dGqN3Ouro6NQqajJKpOzMohVQ4US21ZSpjm6N0R22jo5OreIvHCSSCJSd7L4OkFRNGUs1dGuMHsPW0dGpVfyFo0S2A02EEA2FEBbgdmDBBTILgHtKR4t0AnKllMnV1K0xegtbR0enVvFXTf4kpXQIIR4DlgFG4FMp5QEhxMOl6TOBxcAg4ChQCIz0pPtnbdIdto6OTq3ir5w0V0q5GM0pl903s8y2RFvcpVq6fxbdYevo6NQqZIX9fbUD3WHr6OjUKhz/X+fD1tHR0bnc0FvYOjo6OpcJ/46F3/4eLrnD3ngijakr/0CVkiFt6jLq6ji39O2nM3nq5x3EBGork/dpGsVDXZoA8M3OE/y09zRSws1t6nF3+4YAfLDhEGuOpCKEIMTHwquD2hLhZ/Voh9fVHQh48jEwGCn8dREFX3/rlm6sV5egieMxN21C/idzKPh2njPNd+iteN9wHUiJ4/hxcl5/A0rsWHv1wG/UCEz165H5wGjshw57tMG7S3tCnn0EDAbO/byE3M++d0s3N6hL6CvP4NUijuwZn5H35f802yLDCXvtWYyhISBV8n9cTP7cnzWdpo0Iff5JDD7eOJJSSJ84BVlQ6NEOgCteu4fIPu1QbCX8/uRMcvedLCfTcFQ/Gj8wAL+GUSxu+ZBzVfKo/vG0GH8bqCqqorLvxa/I2nYIgEb3D6DB3b1ACE59vYpjnyytsPwXXp/Guo3bCAkO4pevZ5ZLl1Iy+d2ZrN+8HavVi/8+/zQtm2nXzoYtO5jy7kwUVeWWGwZw//ChAOTm5fP0i5NJSkklJiqStydNIDDA32M9bDx8lqmLdqCqkiHt4xjVo3U5me3HU3hz0Q4cqkqwjxdzHujPyfRcnv1uvVPmbPY5Rvdpy93XtOC3faeYuWoPJ9Jz+frhQbSqE+rRBlObDngPf0xbnHnNYooXul+bpvgueN86EqREKgq2rz5AObzfo67P4y9ijNaGBQsfP2ThOfInPognjI3aYOk/HIQBx+412DctrFDOEN0I68j/UPzTdJSD2qK/po4DMF/ZE6RETU+keMEsUOwYW3TE0v1mRFgMRZ++jJp8wqMNF4Pewv6LUFTJ5OUHmDn0aiL9rdz11QZ6NI6kcZj7j+jKOiFMv6WD276j6fn8tPc0X9/dFbNR8OgP2+jWOIL6wb7c26ERj3ZtBsDcnSeYtekIL/S7onJDDAYCxj5J1lPjUNLSCZs9k+INm3CcPOUUkXn55L07HWv3ru6qYWH43Hoz6XePgJISgl59Ge8+vbEtWYbj+AmyJ75E4LNjq64Mg4GQCY+T+vB4HKkZxHwzg8K1m7EfP+2qr9x8sqZ+gE+vay6oSIXstz+m5OBRhI83Md9+SNGWndiPnybs5bFkTZtF8c69+N3Un8B7byPnwy88mhLZpx1+jaJY0XkswVfF0faNUawb9FI5uaxth0hd/jtdf3rRbX/6+v2kLNNWjg9oUZcOs55kZbdn8G9ehwZ392LtwBdRSxx0/vY5UlbspuBESrm8Bw/qy5233MjESW9VaOP6zds5nZjE4u/nsPfAQSa9NYNvP3kXRVF47e0P+OTd14mKCGPY/U/Sq+vVNG5Yn9lfzaNT+3bcP3wos7+ax5yv5zH2kfsqrQdFVZm8cBszR15LZIAPd320hB4t6tA4Isgpk2crYfKCbXwwog/RQb5knbMB0CA8kHmPX+/Mp98bP9K7peYg4yKDmHZnDybN3+rpNGgIA94jnqRg8jjUrHT8J32E/fdNqGdd16Zj/+/k79wEgKFuI3yfeIn8cSM86hZOn+TUt971MLKwioV/hcAy8F6KvpmCzMvCet+rOA7vRGYklZfrMwzl+F7XLv9gzB37YZs5Hhx2vG5+HFOrTjj2rkdNS6Toh/fwum5U1XVxkdTmFvYlfXFmf3IOdYN9qBPkg9looH/zGNYcTa2W7vGsc7SJDsbbbMRkMBBfN5RVh7Ufvp+Xa44Om12p8v5qbtEcJTEJJSkZHA5sK1bh1dXdKao5OdgPHkI6HOX0hdGI8PICowHh5YWSoa2G7jh1GuXMmXLyFeHVuhmOM0k4zqaAw0HBsjX49OzibkN2DiUHDsMFNigZWZQcPAqALLRhP34aY0SYdmz161C8U/vx2Lb8jk+fblXaEtU/ntPztNZh9u9HMQf44FXGSZ0nd/8pCs+UX8pTKXQt/Wn0sTpn1fFvEkvWzqMothKkopK5OYHoQe0rtKF9uys8tn5Xb9jCjQP6IISgbesW5OefIz0ji30Jh6lXJ4a6sdGYzWYG9unBqvVbNJ31m7lp4LUA3DTwWlat2+yxHvYnZlI3xJ86If6YTUb6t6nPmgT387lkzwl6t6pLdJAvACF+3uXy2XoshToh/sQE+wHQKCKQBuHVm6vD2Lg5aupZ1PRkUByUbFmFOd79uqC4yLkpvFz1XS1dwHJ1T+ybPM9fYohpjJqVisxJB1VBObAFU9P4cnKmDv1wJGxHFuRdkIERTBYQBjBbkOeyAZCZScis5OpUxUWjIKr9udyoscMWQnx5sYWlnSsiyt91gUf6W0k7V1RObm9SNkM/X8ej/9vG0QztsTsuzI+diVnk2Eqw2RU2HE8jNd/m1Jm+/iD9Z65kccJZRndt6tEOY3gYSlqa87uano4xPKxax6BmZHDuu3lE/Pg9Eb/8iCwooGT7jmrputkQEYYjJd353ZGa4XS6NcEUE4mleRzF+w4CUHLsJN49OwPg27c7pqjwKvPwjg7GlpTl/F6UnIV3dHCN7Ige2J4+69+i89fj+P0pbWqGvINnCOvUHHOwH0ZvC5F92uET4zkcUBmp6ZlElamfyIgwUtMzSEvPICoi3G1/Wrp2A83MziE8LASA8LAQsnJyPZaRlldIVKCvK68AX9JybW4ypzLzyLOVcN/s37jjg0Us3HWsXD7L9p5kYJsGNT5GAENIGGpmmWszKwNDcPlzaG7fFf83P8d33OsUznqz2rrG5m1Qc7NRU896tEP4ByPzXNeEzM9C+AeXkzE1a4/j95Vu+2V+NvbNi/F54j18xsyA4kKU4/urOPK/DlVU/3O54dFhCyEWXPBZCNx8/rsHvQeFEDuEEDvmrHM9KlX0jv+FddYiMoAlD/Vm3oju3H5VA576WXOGjUL9GdmxEQ/P28qj/9tG04gAjAaX+Y93a86yh/swqEUs3/1+Co+ICs5UNZepEP5+WLt2IX3oHaQNvhVhteLd79pq6f5VNjiz8LYS/tZLZL35kTNOnfny2wQMu4nouR8gfL2R9vJPCNWxpaardiQv2cHKbs+wdeQ0LZ4NnDuSxJEZC7nm+wl0njue3AOnUB1KzTJ22lPeICFEhXZWVLXVK6PqvBRFkpCUxYx7evHhiD7MWr2PUxmu1qXdobD2YCJ9r6h/cUZU1OqrwDD7jg3kjxtBwTsvYb1tZLV1LZ17Y99cjdkBq1GJlr53U7Lqu/L2WX0wNbuKwhlPUfje42D2wtj6mooz+RtQEdX+XG5UFcOuA/wBzMY1A1V74G1PSmVnwLLNHus8m5F+VlLKtIpT84sIv6BzsGx4o1ujCF5fvp/swhKCfSwMaVOPIW3qAfD+uoNE+pfvWBzYIobHf9rOIx5a2UpaOsaICOd3Q3i4M6xRFV7t41GSU1BLW2tF69ZjvqI1tt9WVEvfaUNqulvr1xQZhpJePRs0BSMRb79MweJVFK7a4NxtP3mG1NHa4hamerH4dKt4RseGI/vS4K5eAGTvPo53TIgzzRodQlFKdk0Ox0nmloP4NojAEuJPSVY+p75dw6lv1wDQYsIwipJrcIxliIoIIyXNFY5JTcsgIiwUu8NBSlq62/7wMK0VHxocRHpGFuFhIaRnZBES5DksERnoQ0quK7abmldAeIB3OZkgXy+8LWa8LWbiG0RwKDmb+mEBAGw4nETzmBBCKwiVVAc1Kx1DaJlrMyQMNad8GOo8ysG9GCJiEH4BVesaDJg7dCX/hYertEPmZSECXNeE8A9B5rtfE4aYhngNeUxL9/HHFNeWYlUFoxE1Jx0K80tt3IGxThOU/RurLPevoBavEFZlSKQ9sBN4Hm1SkzWATUq5Vkq5tqaFtYoO5HR2AWdzCrErKssOJtEjLtJNJuNckbM1tS85ByklQd6aE88q0GKlyXk2Vh1JYWCLWABOZbt+ZGuPpdIwxM+jHfaDBzHWjcUYHQUmE97X9qZ446ZqHYOSmoa5VUvw8gLAEn+VW2dldSk+cAhTvVhMMZoNvv17UrjWc4y1LGEvP439xGnyvv7Rbb8huDT2LARBD9xF/g+/Vqh/4rPlrL52IquvnUjy0h3UG6rFuoOvisORb6M4Lafatvg2cJ3DwCsaYDCbnCNILKWOzDs2lJhBHUj8ufrHWJaeXTuxYOlKpJTs2Z+An58v4WEhtG7elNOJSSQmpWC321myci29unZy6sxfot1I5y9ZQa9unT2W0So2lNOZ+ZzNysfuUFi29xQ9mtd1k+nZoi67TqbhUFRsJQ72ncmgUUSAM33p3hMMuMhwCIBy/CCGqFgM4VFgNGHp1Bv7Tvc6M0TGOLeNDZogTGbkubwqdU2t41GTziCzKr8BnEdNOo4hJAoRFA4GI8ZWnXAc/t1NxjZjLLYZT2Gb8RSOhG0UL/kc5fBOZG4mxtg4LYYNGBq2Qs3wHIL5K1Fr8Lnc8NjCllKqwDtCiB9K/6dWpeOxMIOB565tzej/bUNVJTddUYe4MH9+2K05vNva1WfF4RTm7T6FySDwMhmZcsOViNLHs6fn7yS3yI7JIJhwbWsCrJojf3/tQU5mn8OAIDrQm+f7ehghAqCo5E17n5BpU8FgwLZoCY4TJ/G56QYACucvxBASTNjsjxG+PqBKfG+7lfS7R2D/I4Gi1WsJ/3QWUlGwHz5C4QLNKXp170rgmCcwBAUS/OZkHEeOkfX0s5XakDVlBpEfTdaG9c1fhv3YKfxv1UYa5P/vV4yhwUTP/QCDrw9IScBdN3P25vuxNGmI3w19KTl8nJjvtSFw2dM/xbZhG74DexEw7EbtOFZu4Nz8ZVWel9QVu4ns046+W97BYStm15iPnWmdvnmW3WNnUZSaQ6P7+tPk0evxigii16oppK7cze6nPyHm+o7Uva0b0u5AKbKz/aHpTv2Os8dgCfFD2hX2TPgMe27FoxPGvTyF7bv2kpOTR5/Bd/PIfcNxlHa2DhtyHd07d2D95u0MHDoKb6uVSROfAsBkMjLxqdE8NPYFFEVhyPX9iGukhSPuHz6Up198nZ9+XUZ0ZDjTXnveYz2YjAaeu6Ejoz9fiSolN10VR1xkED9s1YZn3nZ1UxpFBNKlaQxDp/+KEDCkfRPiIrXYrq3EwZajybwwuJNbvqsOnGbKr9vJLiji8S9X0Sw6mI9GVhJGU1Vsn0/Hd/wbYDBSsnYJ6tmTWPpo12bJyoWYO3TH0q0fKA5kSTEF01/1qHseS+delFQnHAIgVUqWfoH1jmfBYMCxey0y4yymq3oD4Pi98nzUpGM4Erbhff9roCqoqadw7FoNgLFZeyz970H4+GMd9gxK6imKv51aPZuqiXqxMbHLAFFRbLBSYSGuA66RUk6srk7ZkMg/xb9nxZl//j2l3SlVd0JeCvQVZ1z8O1aciala6BLg+8LXf9rbfh99V7V9zrDkby4r714jDyKlXAQs+pts0dHR0fnTXI6jP6rLP9/k09HR0fkLuRxHf1QX3WHr6OjUKv7xGOzfiO6wdXR0ahV6SERHR0fnMuFyHK5XXfRFeHV0dGoViqj+588ghAgRQiwXQhwp/V9uPgchRF0hxGohRIIQ4oAQ4skyaf8RQpwVQuwu/QyqqkzdYevo6NQqLuGLM88BK6WUTYCVpd8vxAE8LaVsAXQCHhVCtCyT/o6Usl3pp8r1H3WHraOjU6u4hA77JuD83MVfAIMvFJBSJkspfy/dzgcSgNiLLfDvj2GHR/3tRVSFYv939EJYfKoxEdPfjI/8d0T4/g0vrZh63vlPmwCA+kfCP20ChqaeZ7i8nKjJko5CiAeBsis5zCqdC6k6REopk0FzzEKICE/CQogGwJVA2YnRHxNC3APsQGuJe5zER+901NHRqVXUpElSdqK6ihBCrAAqanV6nuegfD5+wI/AGCnl+ekdPwImoY1EnIQ2qZ7HlR10h62jo1OruLgJfCtGSlnp3MlCiFQhRHRp6zoaSKtEzozmrL+RUv5UJu/UMjKfABXP1FYGPYato6NTq7iECxgsAO4t3b4XmH+hgNBmrpsDJEgpp12QFl3m6xCgylUedIeto6NTq7iEnY5TgL5CiCNA39LvCCFihBDnR3xcAwwHelcwfG+qEGKfEGIv0At4qqoC9ZCIjo5OreJSdatLKTOBPhXsTwIGlW5voMKlgEBKObymZeoOW0dHp1ahzyWio6Ojc5mgzyWio6Ojc5nwV44S+behO2wdHZ1ahVqLgyKX3GFvPJTI1PlbUKXKkI7NGNWrbTmZ7ceSeXPBFhyqSrCPlTmjrwPgmw37+WnrISRwc8dm3N2tNQDTft3GuoTTmI0G6oQG8MrQbgR4e3m0w6tTB4LGPoYwGChYsJj8L791SzfVr0vwi89iadaE3Jmfcu6bec60qJ/nIgsLkaoKikLaiNEAhLz2Iqb62qKtBj8/1HPnSBv+IJXhdXUHAsc8BkYjhQsXce6r8jYEPT8ec9Mm5H08h4JvNRuM9eoS8upLTjljbDT5n3xGwbwf8X9gJNZu14AqUXKyyXntDdRqrAjf7L/3Et7nShRbMfuf+Ij8fSfLyXjXC6fNx09iCvIlf99J9j06A2lXMPl7c8WHj2GNDUMYDZz86FeSvtPWaK73wEDq3K2tA5j4zSpOz1pSYfkbD59l6qIdqKpkSPs4RvVoXU5m+/EU3ly0o/S68GLOA/05mZ7Ls9+td8qczT7H6D5tufuaFvy27xQzV+3hRHouXz88iFZ1Qj3WwQuvT2Pdxm2EBAfxy9czy6VLKZn87kzWb96O1erFf59/mpbN4gDYsGUHU96diaKq3HLDAO4fPhSA3Lx8nn5xMkkpqcRERfL2pAkEBvh7tMPYuA2W/sO1tRR3rcG+cWGFcoaYRlhHvULxj9NRErYhQqPxuuVxV3pwBCVr/odj61IATB36Ye7QF1QVx9Hd2Fd8W2G+59l4JEk7J1IyJD6OUd1blZPZfiKVNxfvxKGoBPt6Mee+vgDk2Up49ZctHE3LRQD/GdKJtvXC+WDFHtYcTEQIQYivF6/e3JmIAB+PdtSUf8e7vH8Pl9RhK6rK5J83MfOBAUQG+nLX9AX0aFmPxpGuSa7ybMVM/nkTH9zXn+hgP7LO2QA4mpLFT1sP8fXjN2E2Gnh0zjK6Na9L/fBAOjWN4YmB7TEZDby7eBufrt7DmEEdKzfEYCB43JOkPz4OJS2diM8/wrZ+E44TrtXP1bx8ct6egXePayrMIv2Rsai5eW77sl6Y5NwOfOJh1IKKF5w9b0PgM0+S+aRmQ/icmRSt3+S2Arual0/uO9Oxdu/qXo+nz5A+4gFnPpHzf6Bo3QYAzn3zPfmffAaA72034z/yHnLffKdyO4CwPu3wbRjNhk5jCIytcy70AAAgAElEQVSPo+XU+9k68IVyck1euJNTHy8i5ZfNtJh6H7F39ibxi+XUHdWfc4fOsmv4m5hD/em68R2Sf9yAb+No6tzdmy0DnkeWOLjquwlkLN9F4YkU9+NRVSYv3MbMkdcSGeDDXR8toUeLOjSOCHLK5NlKmLxgGx+M6EN0kK/zumgQHsi8x6935tPvjR/p3VK7acZFBjHtzh5Mmr+V6jB4UF/uvOVGJk56q8L09Zu3czoxicXfz2HvgYNMemsG337yLoqi8NrbH/DJu68TFRHGsPufpFfXq2ncsD6zv5pHp/btuH/4UGZ/NY85X89j7CP3VW6EEFgGjqDo68nIvCys90/Cceh35IWrjguBpc/tKMf2OnfJzGSKZk10pns/NQPl4A4ADA1aYmoWj+3jCaA4wCcAT2jnZDszR/TWzsnMpfRoXofGEYFOmTxbCZMXbuODe3qXnpMiZ9rUxTvo0iSGt+7ojt2hYLNrgYp7u7bk0Wu1RtrczQeZtWYfL9x4tUdbakrtbV9f4nHY+8+kUzcsgDqhAZhNRvq3bcSaA6fdZJbsOkbv1vWJDvYDIMTPG4Djabm0qReBt8WEyWggvlEUqw5ozq1L0zqYjNqhtKkXQWpOoUc7LC2b40g8i5KUDA4HtuWr8O7exU1Gzc7BnnAIHBcXEfO+tie23ypfXNXcsjmOxCSXDStWaS3jCm2ofA4Sr/ZXoZxNQknRXpqSha5jF1YrVGOR5fAB7Un6YR0AuTuPYgrwwVLGWZ4npGsrUhdqzi9p3joiBrbXEqTE5GcFwORrxZ5zDulQ8W0SS87OI6i2EqSikr0pgYhBHcrluz8xk7oh/tQJ8deuizb1WZNwxk1myZ4T9G5Vl+ggX82W0uuiLFuPpVAnxJ+Y0munUUQgDcIDy8lVRvt2V3hs/a7esIUbB/RBCEHb1i3Izz9HekYW+xIOU69ODHVjozGbzQzs04NV67doOus3c9NA7WW5mwZey6p1mz3aYIhtjJqdisxJB1VBObAFU7P4cnKmjv1xJGxHFuRVkAsYG7ZGZqchczMAMMf3oWTjAs1ZAxRWrHee/YmZ1A0tc06uqOCc7D1J75Zlz4l2DZwrsvP7yTSGxDfWyjYZCfC2AOBnNTv1bSUOxN+wnNclHId9yalRC1sI0RXoCOyXUv5W08LScguJCvR1fo8M9GHfmXQ3mVMZeTgUlftmLqKw2M6dXVtxQ3wT4iKDmbF0BzkFRXiZTWw4eIaWdcqvAP7L9sP0b9vIox3GiDCUVNdbpEpaBpZWLWpwJJKw998EJAU/L6TgF/d1iS3t2qBmZeM4c7ZidcAYfoEN6elYWtbEBg3va3tTuHyl2z7/h+7DZ0A/1IICMh+rciw+1ugQis66wiZFyVlYo0MoSctx7jOH+OPIK0Qq2mVelKTJAJyes4wrvxpHj70fYfTzZu+D74GUnDt4hrgJt2MO9kMpKiHs2nbk7Tlervy0vAuuiwBf9p3JcJM5lVl6Xcz+TbsuujTnhisbu8ks23uSgW0aVHm8F0tqeiZREWEuOyPCSE3PIC09g6iIcLf9+w4cAiAzO4fwMK2ewsNCyMrJ9ViG8A9B5rrOhczLwhDb+AKZYEzN21P05X+xxFQccjO26oRj/yaXTmg0xnrNsfQeCg47JcvnoiaVPxfnScuzERXoClVEBvqwL9E9tOY8J3OWU1ji4M5OzbjhykYkZucT7GvlpZ+3cDg5m5axITw7qD3eFs3dTF++m193n8DPauaTUZW++X3ROETtbWN7bGELIbaV2X4AmAH4Ay8LISqa+/W87INCiB1CiB1zlrkeRyuqxgvvsIqqknA2gxmj+vHh/QOYtWI3p9JzaRQZxMiebXj4k6U8OmcpTaNDMRrcdT9ZuRujwcCgC37IFZVajmq0RM+T9sATpN37EBljnsP31sFY2rVxS/fp15tCD63rymyQNbABAJMJr65dKFq11m13/sdzSB0yDNuyFfjeMqRmeVZmS4VVpsmE9WpL/v5TrG0zms29x9Ni8kiMft4UHEni5IwFxM97nvhvJ5B/4BTSUb5dU9FhiwvKUxRJQlIWM+7pxYcj+jBr9T5OZbhaiXaHwtqDifS9on6Nj7W6VHR+hBDVsv9Pluz2zdJ/OCUrvqv8mjUYMTWLx/GH67cnDAaE1ZeiOS9TsnyuW6y74hIrONYLvitq6TkZ3osP7+nFrDX7OZWRh6JKDiZnMbRDE75/dBBWs4lP1x1w6j3etx3Lxg1hUJsGfLflsOdDvwhkDT6XG1WFRMxlth8E+kopXwH6AXdVpiSlnCWlbC+lbH9ff1d8KjLQh5RcV1w3NbeQ8As6HCIDfenStA7eFjPBvlbiG0VxKDkLgCEdm/HdmMF8Ovp6Any8qBfmisMt2HGE9Qmnef2Onogqfi1KWjrGSNdMiMaIMJSMDA8a7pzvxFOzcyhaswFLq+auRKMB715dsa1Y7dmG9AtsCA+vVudgWaydr8Z++DBqdsUzMtqWr8Taq3uFaXVH9qPTyil0WjmF4tRsrLGuDjlrdAjFKe552jPzMQX4IEpDT9YYl0zM7T1IXaTd220nU7GdTsO3SQwAZ+euZkvfCWwf/Ar2nAIKjyeXs6XcdZFXQHiAdzmZLk1iXNdFgwgOJbts3HA4ieYxIYRWECr5q4iKCCMlzXWdpKZlEBEWSmREGClp6W77w8O0+gwNDiI9Q7t+0zOyCAnyHKKR+VmIQNe5EAEhyPwcNxlDdEO8bnkM7yfexdSyI16DRmAsEzYxxrVDTT4JZcIlal4WjoPbte2k45qz96k8/BMZ4ENKriu8lppbSLi/dzkZ7ZyYXOckJZvIAB8iAny4oq72NNK3VT0SSn/DZRnYtgEr/zhdbv+fpTaHRKpy2AYhRLAQIhQQUsp0ACllAdpKCjWiVZ1wTmfkcTYrH7tDYdme4/RoWc9NpmfL+uw6mYJDUbGVONh3Oo1GpR0d5zuakrPPsWr/SQa201rSGw8l8vmavbw7oq/zscsTJQkHMdWNxRgdBSYT3n17Y6sitngeYbUifLyd215Xt8d+7IQz3atDPI6TZ1DSPN8A7AkHMdUpY8O1vSnasMmjzoV49+2Nbbl7S95YxzU3urVrFxynKv5BnPnsN7b0eY4tfZ4jbckOYm7THHtgfByO/EK3cMh5sjb+QeQN2g04Zmh30pdqHVpFZzMJLR2xYwkPxKdxDLZTWrjHUnpTtcaGEjmoA8k/lz/GVrGhnM7Md10Xe0/Ro3ldN5meLeqy62Sa67o4k0GjCNcNe+neEwz4G8MhAD27dmLB0pVIKdmzPwE/P1/Cw0Jo3bwppxOTSExKwW63s2TlWnp17eTUmb9kBQDzl6ygV7fOHstQzx7HEBKFCAoHg1ELbRze6SZjm/4UtvfHYHt/DI4/tlG8+HOUQy4ZU+vObuEQAOXQTowNtYVOREgUGE1QmF+pHc5zkn1OOyf7TtGjeR33+mheh12nypyTxAwahQcS5u9NVKAPJ9O1G8bW4yk0Ku1LOJXpuomsPXiWhmGeOz8vBhVZ7c/lRlXeLRDYifY0JIUQUVLKlNK5XWv80GcyGnjups6Mnr0UVZXc1KEpcVHB/LBZm8D9ts4taBQZRJemdRj6zs8IobWq46K0GODTX64kt7AYk9HAhMFdCPDRhu5N+WUTJQ6Vhz/Rhi+1qRfBC7dUPLoDAEUl563phL3/BsJgpGDhEhwnTuI75AYACn5eiCEkmIgvZmLw9QFV4nf7LaTePhJDYCChU18FQBiNFC5bSfGW7c6sffr2qkY4RLMhd9r7hL4zFYwGCn/VbPAZrNlQ+ItmQ/inHyPO2zDsVtLuHIEsLER4eeHVIZ6cN9wmACNg9IPa0EJVRUlJJWeq5xEiABkrdhHWpx1dt76HYivmwJOuIW1XfjOeP8bOojg1myOvzaXNx08Q99ww8vadJHGu9hRxfNpPtHp/NJ3XTEUIwZFJc7Fnac6g7ZyxmIP9kA6FhAmf4cgtP3LGZDTw3A0dGf35SlQpuemqOOIig/hhq/a4fNvVTWkUEUiXpjEMnf6rdl201/o1QOu82nI0mRcGd3LLd9WB00z5dTvZBUU8/uUqmkUH89HIymOm416ewvZde8nJyaPP4Lt55L7hOEo7fIcNuY7unTuwfvN2Bg4dhbfVyqSJWv+AyWRk4lOjeWjsCyiKwpDr+xHXSAvN3D98KE+/+Do//bqM6Mhwpr1WxTTKUqVkyedY7xoPwoBj91pk+llM8dqUFY6dKz3rmywYG7WmeNEct92OXWvwuvFBvB+eglQcFM8vP2zRLRujgeeub8/oL1Zpv9WrGmvnZFvpOelYek6axDD0g0UIIRgSr503gPHXtWfi/zZiV1Rig/149Wbt3Lz/225OZuRhEILoIF+ev9HDaK6L5PJzw9VH1DhuCgghfNBWWzhRlaxt/tR/vP4yX1/6T5sAgMH4j1cF+49F/tMmANDtg/LjrC81/5YVZ0o+fPGfNgFD82b/tAkAeA996U9H/59pcEe1f2hvnfz2snqR/aLGYUspC4EqnbWOjo7OpUapxW1s/dV0HR2dWsXl2JlYXXSHraOjU6uoaEhibUF32Do6OrUKvYWto6Ojc5lwOQ7Xqy76mo46Ojq1ikv1pqMQIkQIsVwIcaT0f3AlcidL127cLYTYUVP9sugOW0dHp1bhQFb78yd5DlgppWwCrCz9Xhm9pJTtpJTtL1If0B22jo5OLUPW4O9PchPwRen2F8Dgv1v/b49hK5u2/N1FVElGit8/bQIAjQeW/NMmEJXkYY7uS0jxz9V4G/RvRv0j4Z82AQDLI5OqFvqbcexZ8U+b8JdRk05HIcSDaPMknWeWlHJWNdUjpZTJAFLKZCFERCVyEvhNCCGBj8vkX119J3qno46OTq2iJi3nUudZqYMWQqwAoipIqmKOATeukVImlTrk5UKIg1LKdTXQd6I7bB0dnVrFXzmsT0pZ6eQzQohUIUR0aes4GkirSE5KmVT6P00I8TPamgLrgGrpl0WPYevo6NQqFCmr/fmTLADuLd2+F5h/oYAQwlcI4X9+G21q6v3V1b8Q3WHr6OjUKi7h9KpTgL5CiCNA39LvCCFihBCLS2UigQ1CiD3ANmCRlHKpJ31P6CERHR2dWsWlejVdSpkJ9KlgfxIwqHT7ONC2Jvqe0B22jo5OrUJ/NV1HR0fnMqE2v5quO2wdHZ1ahT5bn46Ojs5lwl8w+uNfyyV32MamV+J14ygQBuzbV2Bf83OFcoY6cXg/OpmiudNQ9m1GBIbiNewJDP7BSKni2Loc+8ZFTnlzl0GYuwxEqgpKwk5Klnzl0Q7/HlcR+/L9CKORzO9+I+2jH8vJxP7nAQJ6tUe1FXP6mXex7T+OOTqMeu+MwRwejFQlmXOXkfHZQje98AcHE/v8KPa1uwslu/KFTo0t47EOHQ0GA/aNSylZNs8t3dS2E5Yb7gWpgqpQPO9jlGMHtERvX6zDx2CIaQBSUvTlO6gntDf3zD1vxNLzRq0u9m+j+Kc5eMKvx1XEvvQAGA1kfb+c9I/+V04m5uUH8e8Vj2orJvGZ97AdOIbwMtP4+ykILzPCaCR3yUZS35kLQOTYuwjoezVIiSMjlzPPvIsjrfzK2c5jbdMB7+GPgcFAyZrFFC/81j09vgvet44EKZGKgu2rD1AO7/eo6/P4ixijtcV8hY8fsvAc+RMfxBPGxm2w9B8OBgOOXWuwb1xYoZwhphHWUa9Q/ON0lIRtiNBovG553JUeHEHJmv/h2KoNCDB16Ie5Q19QVRxHd2Nf8W2F+QK88Po01m3cRkhwEL98XX7tRSklk9+dyfrN27Favfjv80/TslkcABu27GDKuzNRVJVbbhjA/cOHApCbl8/TL04mKSWVmKhI3p40gcCAyldNB9h44ARTf1iNKiVDurRmVP+r3dI/X76dxdu1a05RVE6kZLF66mgCfb0Z+MIn+FotGAwCk8HA3OfuBuDZ2Qs5maatdp9fWIy/jxfzJt7j0Y6aoodE/iqEAa/BD2Cb/QoyNxPvx6bi+GM7Mi2xnJxl4HCUw7td+1SVkl+/QE06DhYrPk+8hePIHmRaIsZGrTG27EDhO0+B4kD4Bnq2w2CgzqSHOHbXS9hTMmm64G1yV2yj+MgZp4h/r3i8GsaQ0OMhfK5sRp3XRnNk8DikopD02qfY9h/H4OtN01+nkb9ht1PXHB2Gf9d2lCRWMQZeGLDe8SiF701EZmfgM+F9HHu3oCa7Vjl3HNyNY4/2ar8htiHWByZS+J8HALAOfRjlwE6KZv1XWwHboi1IbGzaBlPbzhS8NhocdoR/1XUR++rDnLj7RewpmcQtmEbe8q0UHy1TFz3jsTSM4VBPrS5i/zuao4OfQRbbOX7n86iFRWAyEve/N8hfs5PCXYdIn/UTqdO+ASB0xA1EPnk7Z5//sNK68B7xJAWTx6FmpeM/6SPsv29CPXvKVRf7fyd/p7YSuKFuI3yfeIn8cSM86hZOd73ybb3rYWRhFa/lC4Fl4AiKvp6MzMvCev8kHId+R2acLS/X53aUY3udu2RmMkWzJjrTvZ+agXJQm5jN0KAlpmbx2D6eAIoDfDyvFD54UF/uvOVGJk56q8L09Zu3czoxicXfz2HvgYNMemsG337yLoqi8NrbH/DJu68TFRHGsPufpFfXq2ncsD6zv5pHp/btuH/4UGZ/NY85X89j7CP3VWqDoqpM/n4lM5+4lcggf+564xt6tImjcXSoU2ZE3w6M6NsBgLV7j/H1qp0E+no70z8ZcxvBfj5u+U69/wbn9ts/rsHP28tjXVwMtbnT8ZKOwzbUjUPNTEZmpYLiwLFnA6aW5VdNNl8zCGX/ZuS5XOc+mZ+tOWuAkiLUtEQMgdrFY+rcX2upK9oK17Igt1yeZfFp14Tik8mUnElF2h1kL1xPYF/31kNg36vJ+lFbFbxw1yGMAb6YIoJxpGVj26/ZoRbYKD6aiDnSdRHHvnQfSZM/hyoeywwNmqGmJSMzUrS62L4WU5vO7kLFRa5ti9WVp9UHY5MrsG8sHc6pOMCmOSNzj+u1lrrDXlpvVddFySlXXeQsXEdAP/e6COjXiZyfVrnqwt8XU7g2E6RaqNkoTCaEycT5RZ3VczbXsfp4eawPY+PmqKlnUdOTQXFQsmUV5vguldaF8HLVRbV0AcvVPbFv8jx/iSG2MWp2KjInHVQF5cAWTM3iy8mZOvbHkbAdWZBX8fE0bI3MTkPmZgBgju9DycYFzuuTwor1ztO+3RUeW7+rN2zhxgF9EELQtnUL8vPPkZ6Rxb6Ew9SrE0Pd2GjMZjMD+/Rg1Xrthr96/WZuGqi9tHfTwGtZtW6zRxv2n0yhbngQdcKCMJuM9I9vxpo9RyuVX7LjIAPaN/eYZ1mklPy281CNdKqd96Wb/OmS47GFLYS4GkiQUuYJIbzRpv+7CvgDeF1K6dkbXJhfYCgyJ9P5XeZmYqjXxF0mIARTq6uxzXoZr1vjKs4nOBxDbEOU04cBMITFYGzYAkv/O8Fhp3jRF6iJlV9c5qhQ7MkZzu/25Ax8rmxWXiYp3SWTkok5MhRH6eMcgKVOBN6tGlG4+xAAAdd2xJ6SSVHCySpqAgzBoajZrvzVnAyMDcuvXG1q1wXL4JEY/IMonPFS6fFGIc/lYr336dJ6OErxvI+gpBhDRCzGuFZ43XQv0l5C8Y+zUU8drrwuIkOxJ5Wti0x82jUtJ1NSRqYkJRNzVCiO9GwwGGjy6ztY6keT+dUibLtdZUU+M5zgm3uh5hdy7I6JlddFSBhqpuuJRM3KwNS4RXlb23fFOux+REAQBW9OrLausXkb1Nxs1NQLWsoXIPxDkLllrs+8LAyxjS+QCcbUvD1FX/4XS0zF4RVjq0449m9y6YRGY6zXHEvvoeCwU7J8rqvxcRGkpmcSFRHm/B4ZEUZqegZp6RlERYS77d93QLs2M7NzCA8LASA8LISsHM8/3bScc0QFu24akcH+7DuZXKGsrcTOpj9OMmFYb+c+IWD09B8RwC3d2nJr1zZuOr8fPUtogC/1I6qcArrG1OaQSFUt7E+BwtLt94BA4I3SfZ9VpiSEeFAIsUMIsePT3VUsrn5B3XrdMIriJV9pcduKsFix3v0sxQs+heLSVpzBCN5+2D54juJFX2C96+kqDquCle0vbAFWIFJWxuBjpcHM5zj76mzUczaE1ULkY7eRPG1uFWXXwAbAsXsThf95ANtHr+B1Y2msz2DEUDeOkrW/Uvj6Y1BShKX/MGea8PGn8I0xFP80G+8HKneUmhl/si5UlSODniSh80h82jbFq2k9p0jqW19xsMsosuevIeze6z0ZUbUNgH3HBvLHjaDgnZew3jay2rqWzr2xb77Y2QEvyKv/cEpWfFf5E4PBiKlZPI4/tjp3CYMBYfWlaM7LlCyf6xbrviiLKihbCFGhSRWd3mqVUYHTqyyrdXuP0a5RjFs45POn7+C7CcP54LFbmLd2NzuPuIc9l9awRV4TpJTV/lxuVOWwDVLK0uc42kspx0gpN0gpXwEaVaYkpZwlpWwvpWw/ql1D1/7cTESQK3wgAkORee4dUYY6jbHeMRaf8TMxXdEZr8EPYjwfNjEYsQ4fh2P3OpQDW93yVfZrj35q4lHtx+RbeZzQnpKBOdrVQjFHh2FPdbfDnpyJOcbVWjFHhWI/32lmMtJg5nNk/7KW3KXao6VX/WgsdSNpvuQ9Wm74BHN0GM0WvYspPKhCG9TsDAzBrvwNQWHInMo75ZSj+zGERyN8A5A5GcicDNSTWuvJ8ft6jPW0pxGZk4Fj90atjJOHQaoIv8rj2PaUDMwxZeuizHE6ZTKxlJGxRIWWqy81r4BzW/bh36N8CCFn/loCB5QPUzh1s9IxhLpmljSEhKHmZFQqrxzciyEiBuEXULWuwYC5Q1dKtqyuNL/zyPwsRGCZ6zMgBJmf4yZjiG6I1y2P4f3Eu5hadsRr0AiMZcImxrh2qMknoUy4RM3LwnFwu7addFy7Pn08d/h5IioijJQ01zGmpmUQERZKZEQYKWnpbvvDw7TjCQ0OIj1DO2fpGVmEBHnu24gM8ielTId5anY+4YEVT1O8dOchBnRwd74RQZpsiL8PvdrGsb9M69yhqKzcfYT+8eWfKP8KFGS1P5cbVTns/UKI802ZPUKI9gBCiKaAvaaFqYlHMYRGI4IjwGjC1LYrSsJ2N5nCN0ZT+MbDFL7xMI59myn+ZRbKH9sA8Lr1UdS0s9jXu/fcOw5sxdj4CgBEWLTWCVdJfBGgcM8RvBrGYKkbiTCbCL6hG3nLt7rJ5K3YRsgtvQDwubIZSn6hMxxSb+rjFB9NJH22a66WokOnOBB/D390fYA/uj6APTmDQ9eNwZHu/oN31sWpQ5rTCY3U6qJDDxx73ecOF+HRzm1D3TgwmZAFeci8bNSsdERkHQCMza90dlY6dm/C2Ex7E1ZExILR7NYXUFFdWBrEYK6j1UXQDd3JW77NvS6WbyXo5t7udZGejTEkAEOAr1aWlwX/a9pRfExrSVkauGwPuPZqio5d0LFcBuX4QQxRsRjCo8BowtKpN/ad7jFWQ2SMc9vYoAnCZEaey6tS19Q6HjXpDDKr8hvAedSzxzGERCGCwsFg1EIbh3e6ydimP4Xt/THY3h+D449tFC/+HOWQS8bUurNbOARAObQTY8OWWj2FaHZSWPnooaro2bUTC5auRErJnv0J+Pn5Eh4WQuvmTTmdmERiUgp2u50lK9fSq2snp878Jdqc1/OXrKBXt86eiqBV/ShOp+VwNiMXu0Nh2c5D9GjTuJxcvq2YnUcS6dXGFb60FdspKCpxbm9OOElcmRv+1oOnaBgZQmTwxd+0PHEJ5xK55FQ1SuR+4D0hxAtABrBZCHEGOFOaVjNUleL5s/G+7yVtKNv2laipZzBd3Q8Ax9bfKlU1NGiOOb4nSvJJvJ98G4CSpd+gHPodx45VeN36KN5PvQuKg+J573u2Q1FJfOljGn35H4TRQNa8FRQdOUPoXQMAyPxmKXmrduDfK54W6z4uHdan5enbvgUht/TGlnCSZovfBSDpza/IX72z0uIqq4ui7z/E54n/anWx6TfU5FOYuw0CwL5+MeYru2LqdK3WWWUvoeiTyU714u8/xHvUs2A0o2YkU/TlNE1v029Y7xmLz4szQXFQ9EXFIw3K1kXSSzNp9OUrYDSQPW8FxUdOE1JaF1nfLCV/9Q78e7Wn2dpZ2rC+ce8BYI4Ioe7bY8BgQBgM5CzaQP4q7QYcPX4EXo1ikaqK/Ww6ic9/4LEubJ9Px3f8G2AwUrJ2CerZk1j6aCMKSlYuxNyhO5Zu/UBxIEuKKZj+qkfd81g696KkuuEQqVKy5HOsd40HYcCxey0y/SymeG26B8fOlZ71TRaMjVpTvMh9GKVj1xq8bnwQ74enIBUHxfPLD9Ury7iXp7B9115ycvLoM/huHrlvOA6H9qA7bMh1dO/cgfWbtzNw6Ci8rVYmTXxKK95kZOJTo3lo7AsoisKQ6/sR16g+APcPH8rTL77OT78uIzoynGmveZ7O2WQ08Nyw3oye8SOqqnJT59bExYTxw7o9ANzWXWsUrNp9hM4t6uPtZXbqZuYXMPbjBdqxqyoD2zfnmlauJ+2lf1Nn43kux1BHdRHVObjS6QEboTn4RCllanULODf+5n+89o5+56ha6BLwb1hx5sRSc9VCl4B61/zzK9+Ym0VXLXQJ0FecceHd58GLjLq76FWnb7V9zurE5X+6vEtJtcZhSynzgT1/sy06Ojo6f5rLcbheddFfTdfR0alV6K+m6+jo6FwmXI6didVFd9g6Ojq1Ct1h6+jo6Fwm1OZRIvqajjo6OrWKSzUOWwgRIoRYLoQ4Uvq/3Hv2QohmQojdZT55QogxpWn/EUL8X3vnHR5Vlf7xzzslPSGkJ4BAaBZ6ExRBQBTEFVER2RVRYHXtIq4Fdd1Ff44UePUAABaPSURBVLa1CysiqNhXVFaUJoI0BWmCgPQmIT0BkkDItPP74w6ZDJOGkJkhns/z3Ce3vPee7xwu7z33vee852CFY1fWVKZ22BqNpl7hx+RPjwCLlFKtgEXubW8tSm1XSnVUSnUEumCk9aiYU/qVE8eVUnNPPv9ktMPWaDT1Cqdy1Xo5TYYAM9zrM4BrarDvD+xWSu2vwa5K6jyGLdERNRvVMWGhVefo8CeOvOM1G9UxJWWRgZYAgPW8GvJ0+wFT69Y1G/mBYBi0YulwWaAlnDH8GMNOVkplucvMEpGkGuxvBE6eueJuEbkZWAuMV0od8j3Ng25hazSaesWpxLArZhZ1L145c0XkOxHZXMky5FQ0iUgIcDUws8LuN4EWQEcgC3ippuvoXiIajaZecSqxaaXUVGBqNcerfPUQkRwRSXW3rlOB6qaZGgSsr5jWo+K6iLwNfFOTXt3C1mg09QqXUrVeTpPZwCj3+ijgq2psR3BSOMTt5E8wFNhcU4HaYWs0mnqFH3uJPAcMEJGdwAD3NiKSJiLlPT5EJMJ9/MuTzn9BRDaJyC9AX2BcTQXqkIhGo6lXnIHeH7VCKVWA0fPj5P2ZwJUVto8B8ZXYjTzVMrXD1mg09YozEOoIWrTD1mg09QqdXlWj0WjOEnQLW6PRaM4SdAv7DGJOb0/IFSPdc+Ytwf7j15XamVLTCbv1n5R9+QZO94zTlu4DsXa6FJTClZdB2eyp4LRj7XM9ltadjRFOx4oom/0WqqTyyW9PEHlJF5Ieux0xmzg8cwGFU2f62CQ9fjtRfbrhKi0j65GXKft1NyHNG5H2qidlgLVJKvmvfcChGZ4ePXGjryXpkbHsvPBGnIeqngzY0rE7EaPvBpOZskVzKJv1sdfxkEsuI3ToCGOjtJRjU1/BuX83prQmRD7wpKdOk1Mp/fRdyuZ8DkDooKGEDhoKLif2daso/eCtausCoPlTo4nt3xlXqY1d97/B0U17fWxCmyTReso4LLHRHN20h533vI6yO2h4RTfOeWgEuFwop5O9/3iX4tXbAIjt25HmE0eD2UTux4s4OGmWz3Whbu4L83ndCel9LZKQxvF3nsSV5fubTuaHnZm8MGctLqUY2qUlo3tf4GOzZm8O/567DofTRcPIUKaPGQBAUamNif9bxa7cIwjwz6E96HBOIpO/28iSbRmICHGRoUy8tidJMVWPAP5hy15emPm9oeGitoy+4kKv4+8tXMPcNVsBcDpd7M0u5PsX7qBBZDiDHn+byLAQTCbBYjLx8SM3AfDQtK/Z555EuvhYGdERoXw24eYqNTz+zMss+2E1cQ1j+d+HvnNQKqV49tUpLF+5hrCwUP7vsfGc38aYiHfFqrU89+oUnC4X1/1pIGNH3gDAkaJixj/xLJnZOaSlJPPSU4/SIObMT8TrVM4zfs1gwb8OW4SQQaM4/tFzqKJCwsZMxLFjHSo/09eu/3Cce37x7IpuiLX75ZROeRgcdkKvvQfLBT1w/LIc+8o52JcazsrS7XKslwzFNu/dqnWYTCQ/eScHbn0Me3Y+zb54lZJFq7DtPlBuEtmnKyHNGrFnwFjCOrQh5V93s3/YOGx7D7JvyD3l12m5/H2KF1aYpTslgYiLO2E/WF0feuPciL/eR8nEB3EV5BH9/BTsa37AleFJM+DMzaLkiftQR0uwdOpOxN/GU/zonbgyD1D84Njy6zSY+jn21cuN8tt2xNq9F0UPjAGHHYmJrV4HENuvM2Hpqfx80d1EdW5F+nO3sWnwoz52TR8fSebUbyj46gfSn7+NpBH9yXl/AUeWb2LjAsN5RpzXlNZTx7PhknvBZCL9mb+yZfhEbFkFtJ/3PIXfrqF0x0kzqNfRfeHKzeD4zNcIHTy6xjoAcLpcPPv1Gqbc0o/kmAj+MmU+fc5tTIskzzD6olIbz369msk39yM1NpLCEk+6gRfmruWiVmm8OKI3doeTUrvhOEb1Op+7LjMmrf145TamLtnE41d7O2EvDf9dxJR7ryc5Npq/PP8Rfdq3pEWqp5PBLQO6ccuAbgAs/WU3Hy5eR4PI8PLjb98/jIZR3g+EF8b+qXz9pS+WEBUeWm1dXHPlAP583dVMeKrySZyXr1zDbxmZzP3vdH7Zso2nXpzEJ2+/itPp5OmXJvP2q8+QkpTA8LH30bfXhbRo3pRpH3xGj64dGTvyBqZ98BnTP/yMB+4cU62O34NOr3qmCktrgaswB3U4D1xOnFtWYWndxcfO0u1yHFvXoI6e1Do1mcESAmICawiqxD3s3lZabiLWUKjhlSisfWts+zOxH8gGu4OiOcuIuqynl01U/x4cmWXMkn1843ZM0ZGYE72zJ0b07IDtt2wcmR7nnDThNvL+/Q7UcNOYW56LK/sgrpwscDiwr1hMSLeLvWyc27egjpYY6zt+xRSf6HMdS7vOuHIO4sozBk2FXjGE47M+BocdAFVU/ZsGQNzAbuTNXApAyfqdWGIisSb5OvoGvdpS8I3xcMr9bAlxg7oD4DrmcVqmiNDy3x7VqSWl+7Ip+y0HZXeQ/9UK4q7o5nPdurovVEEmqjCrxt9/gs0ZBTSJj6ZxXDRWi5kr2jVlydYDXjbzftlHv/ObkBpr5GSJiwoDoOS4nfX7chnapQUAVouZmPAQox7CPBMfl9ocCFXP+7p5XzZNEmNpnBBraOjShiUbd1VpP2/ttlOagVwpxbe1mLW8a8d21bZ+v1+xiqsH9kdE6ND2PIqLS8jLL2TT1h2c0ziNJo1SsVqtDOrfh8XLVxnnLF/JkEHGwMEhgy5j8bKVVV7/dPBXetVAUG0LW0TuBWYppQ5UZ1dbJLohqsiTiEkVF2JKa+FjY2nTleMfPkNIWnoF20PYV84l4t7XwG7DuXcTzj2egUHWS4dhad8Ljh+j9MNnqtVhTY7HkZ1fvu3Izie8Q5uTbBJwZOd5bHLysSYn4Mzz5GaJGdyHojlLyrej+l2II6eAsm01v3qb4hJx5Xuu7yrMw9zq/CrtQ/oPxv7zat/9F/fDtmKx57qpTbCc147wEWNQdhulM97EuXt7tVpCUuIoy/TUR1lWASGp8dhzPc7eEheN48hRcBp9XG1ZBYSmxJUfjxvUnXMm3IQ1PoatI436D02Jw3bQc11bViFRnVr5lF+X98WpkFtUSkoDT8s0uUEEmzIKvGz2FxThcLoYM30hx2wO/tyjDX/qlE7GoWIaRobxj1mr2JF1iPMbxfHQlV0JDzH+i72xcAPfbNhLVJiVt0dXnWgp93AJKQ09jjK5YTSb9lX+0Cm12fnx1308Orxf+T4RuOONLxDguks6cH2v9l7nrN91kPiYSJom+aRuPiVy8gpISUrw6ExKICcvn9y8fFKSEr32b9pi3H8Fhw6TmGDcM4kJcRQePnJaGqrij9zCfgr4SUSWi8idIuLbxKuEiglV3lmzs+KBGs8NGXATtsWf+rZQwyKwtOnMsUnjOPbaPWANxdzW0yK1L5lJ6ev34dj8I9auA2oS6Lvv5PIqNalgY7UQ1f9CiuetMMzDQom/40byX/ug+rKruX5VrXJL246E9r/SNxZtsWDtdjG2H5d4Lms2I5HRFD96J6XvTyFy/D9rllKr+vC1qVgfhfNWs+GSe9k++gUjnl3FOZX+xjq8L06Fyj5WnazM6VJszSxk0si+/Ofmvkxdspn9+UU4XYptWYXc0K0V/73rSsKsFt5ZtqX8vHsGdGTB34dyZftmfLpqx2lpOMGyX3bTMT3NKxzy3vgRfProSCbffR2fLd3Aup3e4af5p9gir1JnJf+OIvJ7/3nPKH4cmu53anLYe4DGGI67C/CriMwXkVEiUuX7klJqqlKqq1Kq6+hunhaVKipEYjytMomOQxV7ZxM0pTUndOjdhN/9CpbzuhM66BbMrbtgbt4W1+E8OFZsvDZvW4u5sW9rzbHlRyzn+r52V8SenY8lxdM6sKQkYM8trMTG83yyJCfgyPW0tqJ6d6Vsy26cBUYrNOScVKyNk2k+ezItFr+LJSWBZrNex5xQeUvGVZCHKcFzfVNcIqow38fO3DSdiDv+Tslzj6FKvEMB1k4X4tyzA3XEU4eugjzsPxnxbOeubaBcSIxvKtOUWwbSYeGLdFj4IracQkLTPPURmhqPLdu7PhwFRVgaRILZuGVCUuOx5fhmgixa9SthzZKxxEUbLfVGnuuGpMZhy/FNdeuP+6I2JMdEkH3kWPl2zpFjJEaH+9hc1CqN8BALDSPD6NIsie3Zh0iOiSApJoJ2TYzfO+CCc9ia5ftbB3VoxqJff6taQ2w02YeKPRoOFZPYIKpS2/nrtjOwm7fzTYo1bOOiI+jboSWbK7TOHU4Xizbs5Iou3m+Tv4eUpASycz33a05uPkkJ8SQnJZCdm+e1PzHBiL/HN4wlL9+ok7z8QuJi6ybFrh+Hpvudmhy2Ukq5lFLfKqXGAGnAf4CBGM78lHBl7sEUl4LEJoLJjPmCHjh2rPeyKZ30AKWTxlE6aRyOraspm/cezh3rUEcKMDdqacQqAVPzC3DlHwRAGiaXn29u1RlXQfVxy+ObdhDSLA1r42SwWogZ3JuSRau8bEoW/0SDocao07AObXCVHPUOh1zVh6JvlpZvl+3Yx66ef2Z3v1vZ3e9WHNn57Bt6L878ytPbOndtx5TaGFNSitFS7tUP29ofvWwkIYnIvz/F0defwZWV4XONkF79sa1Y5LXPtnoFlnadjDpKbYxYrKgi31fP7Pfms3HAg2wc8CCF81aTOKwPAFGdW+EoPuYVDjnBkR82E3+VEetPuuFSDs03QjRhzVLKbSLbNUesFhyFxZRs2EV481RCmyQhVgsJQ3pRuGCtz3Xr6r44VS5oFM9vBcUcPFSC3eFkwab99Dm3sZfNpec25uf9uTicLkptDjZl5JOe2ICE6HBSGkSwL894qP60J5v0RMMh7S/wPGiXbjtI84SYqjU0TeG33MMczD9iaFi3nT7tW/jYFZeWsW5nBn3bt/TUUZmdo8dt5esrt+6jZYUH8U/b9tM8OY7khqffM+PSXj2YPX8RSik2bt5KVFQkiQlxtD23Nb9lZJKRmY3dbmfeoqX07dWj/Jyv5hm5v7+a9x19L+lZXRG/Gz9OYOB3auol4vUyo5SyY2Somi0i4ZWfUg3KhW3+DMJGPAQmE44NS1H5B7F0NmJwjvWLqzzVlbkbx9bVhI99GlxOXDn7cfz8PQAh/YZjik81unUdya++hwiA00XOxDdpMv1pMJs48vm32Hb9RuyNxvD/w5/O5eiSNUT16Ub6d9NxlZaR/egrnkoJCyXyok5kP/HGKVeB5wc5OTbtNaKe+DeYTNgWz8N1YB8hl18NgO3b2YQPG4VExxDxV3dOGKeT4odvN9ZDQrF06MLRt7xT6NoWzyXizoeJeeVdlMPO0TeerVHKoUXrie3fmc4rJ+MsLWPXuMnlx8778DF2jf8P9pxD7H/6Q1pPGcc5D4/g6Oa95HxiPCziB/cgcdilKLsD13EbO/72sluviz0TpnH+J08gZhM5ny6mdEcln0Pq6L4wt+lKyBU3IxHRhA1/EGfOfso+eaHKa1nMJh65qit3zFiMy6UY0rkFLZNjmbnaCGEM696a9KQGXNQqjRsmz0FEGNqlJS2TjQ+0Dw/uyoTPf8DudNGoYRQTrzUc1evfbmBffhEmEVJjI3ns6u7VaxjejzsmfYHL5WJIz7a0TEtg5rKNhobeRm+TxRt20vO8poSHej5oFhQf5YG3Zht15nIxqOu5XHxB8/Lj82vxsfEEf3/yOdb8/AuHDxfR/5qbuHPMSBwOBwDDhw6md89uLF+5hkE3jCY8LIynJhj3qMViZsK4O7j9gcdxOp0MvepyWqY3BWDsyBsY/8QzfPnNAlKTE3n56cdqpeVUqc8xbKnux4lIa6VU1QG3WnD06ZsCXnsH3g+OGWeS2x2r2aiO2fqjTw6agNDhrlN/3p9pgmXGGeJTarapY4JlxhlrQvppR7zjolvV2ucUFu/0c4T99Ki2hX26zlqj0Wj8TX1uYeuh6RqNpl5xNvavri3aYWs0mnqFbmFrNBrNWcLZ2PujtmiHrdFo6hVn44CY2qIdtkajqVfU55CInoRXo9HUK/w10lFEhonIFhFxiUjXauwGish2EdklIo9U2B8nIgtFZKf7b40JXrTD1mg09QqlVK2X02QzcC2wrCoDETEDk4FBwPnACBE5keXtEWCRUqoVsMi9XS3aYWs0mnqFv5I/KaW2KqWqT4UJ3YFdSqk9Sikb8CkwxH1sCDDDvT4DuKY2hQb9AtymNQSPjmDQECw6gkFDsOgIBg2/RzOwtsJyyr8BWAJ0reLY9cC0CtsjgUnu9cMn2R6qqayzpYV9W6AFEBwaIDh0BIMGCA4dwaABgkNHMGg4JVSFzKLuZWrF4yLynYhsrmQZUtU1T6LSRMq/V6/uJaLRaDRVoJQ63SQrGUCTCtuNgRNz3+WISKpSKktEUoEa5hXUMWyNRqOpS9YArUSkuYiEADdiZDzF/XeUe30U8FUl53txtjjsqTWb1DnBoAGCQ0cwaIDg0BEMGiA4dASDBr8hIkNFJAPoCcwRkQXu/WkiMhdAKeUA7gYWAFuBz5RSJ6Yieg4YICI7gQHu7erLdAe7NRqNRhPknC0tbI1Go/nDox22RqPRnCUEtcOuakinnzW8IyK5IrI5EOW7NTQRke9FZKt7KOx9AdIRJiKrRWSjW8e/AqHDrcUsIj+LyDcB1LBPRDaJyAYR8Z2s0j8aYkXkcxHZ5r4/6maixOo1tHHXwYmlSETu97eOPwJBG8N2D+ncgRGMz8D42jpCKfWrn3X0BkqA95VSbf1ZdgUNqUCqUmq9e7b6dcA1AagLASKVUiUiYgVWAPcppVbVcGpdaHkA6ArEKKWu8nf5bg37MAZM+E537z8NM4DlSqlp7l4IEUop3xmU/afHDBwELlRK7Q+UjvpKMLewqxvS6TeUUsuAgE4KqZTKUkqtd68XY3xtbhQAHUopVeLetLoXvz/xRaQxMBiY5u+ygwkRiQF6A9MBlFK2QDprN/2B3dpZ1w3B7LAbARWn2M4gAE4q2BCRZkAn4KcAlW8WkQ0YnfwXKqUCoeNV4CEg0JnqFfCtiKwTkUCM8ksH8oB33eGhaSISGQAdFbkR+CTAGuotweywz+iQzvqAiEQBXwD3K6WKAqFBKeVUSnXEGLHVXUT8GiYSkauAXKXUOn+WWwUXK6U6Y2Riu8sdPvMnFqAz8KZSqhNwlFpkfKsr3CGZq4GZgdJQ3wlmh13dkM4/HO6Y8RfAR0qpLwOtx/3qvQQY6OeiLwaudsePPwX6iciHftYAgFIq0/03F5iFEcbzJxlARoW3nM8xHHigGASsV0rlBFBDvSaYHXZ1Qzr/ULg/9k0HtiqlXg6gjkQRiXWvhwOXAdv8qUEp9ahSqrFSqhnGPbFYKXWTPzUAiEik+wMw7jDE5Rj5kf2GUiobOCAibdy7+gN+/RB9EiPQ4ZA6JWiTPymlHCJyYkinGXinwpBOvyEinwCXAgnuYahPKqWm+1nGxRhpGTe548cAE5RSc/2sIxWY4e4JYMIYZhuwbnUBJhmYZTxLsQAfK6XmB0DHPcBH7kbNHuDWAGhARCIwenTdHojy/ygEbbc+jUaj0XgTzCERjUaj0VRAO2yNRqM5S9AOW6PRaM4StMPWaDSaswTtsDUajeYsQTtsjUajOUvQDluj0WjOEv4f0NfN/rIssE4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Correlation matrix\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "Xn = (X_train - X_train.mean(axis=0))/X_train.std(axis=0)\n",
    "n, p = Xn.shape\n",
    "C = Xn.T@Xn/n\n",
    "# or C = np.corrcoef(X_train.T)\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "sns.heatmap(C, annot=True, fmt=\".3f\", vmin=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Compute the linear fit to the prostate cancer data, the standart error assosiated with each coefficient and the corresponding Z score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.concatenate( (np.ones((X_train.shape[0],1)), X_train), axis=1)\n",
    "b_ls = np.linalg.solve(X.T@X, X.T@Y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, p = X.shape\n",
    "y_hat = X@b_ls\n",
    "sigma_square = np.dot(Y_train - y_hat, Y_train - y_hat)/(n-p)    \n",
    "vector_v = np.diag(np.linalg.inv(X.T@X))\n",
    "    \n",
    "std_error_coef = np.sqrt(sigma_square * vector_v)\n",
    "z_score = b_ls/(std_error_coef)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------------------------\n",
      "Term       Coefficient    Std. error    Z score   \n",
      "--------------------------------------------------\n",
      "Intercept         2.46        0.09       27.60\n",
      "lcavol            0.68        0.13        5.37\n",
      "lweight           0.26        0.10        2.75\n",
      "age              -0.14        0.10       -1.40\n",
      "lbph              0.21        0.10        2.06\n",
      "svi               0.31        0.12        2.47\n",
      "lcp              -0.29        0.15       -1.87\n",
      "gleason          -0.02        0.15       -0.15\n",
      "pgg45             0.27        0.15        1.74\n",
      "--------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "dash = '-' * 50\n",
    "print(dash)\n",
    "print(\"{:<11s}{:<15s}{:<14s}{:<10s}\".format(\"Term\",  \"Coefficient\", \"Std. error\",\"Z score\"))\n",
    "print(dash)\n",
    "for k in range(variables.shape[0]+1):\n",
    "    if k==0:\n",
    "        print(\"{:<10s}{:>12.2f}{:>12.2f}{:>12.2f}\".format(\"Intercept\", b_ls[k], std_error_coef[k], z_score[k]))\n",
    "    else:\n",
    "        print(\"{:<10s}{:>12.2f}{:>12.2f}{:>12.2f}\".format(variables[k-1], b_ls[k], std_error_coef[k], z_score[k]))\n",
    "print(dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = \"https://i.stack.imgur.com/epSz4.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Best subset "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation:  [ 2.477  0.740  0.316  0.000  0.000  0.000  0.000  0.000  0.000]\n"
     ]
    }
   ],
   "source": [
    "X = np.concatenate( (np.ones((X_train.shape[0],1)), X_train[:,0:2]), axis=1)\n",
    "b_bs = np.zeros(p)\n",
    "b_bs[0:3] = np.linalg.solve(X.T@X, X.T@Y_train)\n",
    "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})\n",
    "print(\"Estimation: \", b_bs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ridge\n",
    "\n",
    "It differs from the regession by the use of a $\\ell_2$-norm regularization on the parameters\n",
    "$$\n",
    "\\min_{\\beta_0, \\, \\beta} \\sum_{i=1}^n \\left(y_i - \\beta_0 - x_i^\\top \\beta \\right)^2 + \\lambda \\sum_{j=1}^p \\beta_j^2\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation:  [ 2.452  0.420  0.238 -0.047  0.162  0.227  0.001  0.041  0.132]\n",
      "2.4641208076992855\n"
     ]
    }
   ],
   "source": [
    "ybar = Y_train.mean()\n",
    "stdy = Y_train.std()\n",
    "xbar = X_train.mean(0)\n",
    "stdx = X_train.std(0)\n",
    "Xc = X_train - xbar\n",
    "Yc = Y_train - ybar\n",
    "lam = 24.25\n",
    "b_R = np.linalg.solve(Xc.T@Xc + lam*np.eye(p-1), Xc.T@Yc)\n",
    "b0_R = ybar \n",
    "# b0_R = ybar - xbar.T@b_R # should be this\n",
    "beta_R = np.concatenate((np.array([b0_R]), b_R), axis=0)\n",
    "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})\n",
    "print(\"Estimation: \", beta_R)\n",
    "b0_R = ybar - xbar.T@b_R\n",
    "print(b0_R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation:  [ 2.464  0.420  0.238 -0.047  0.162  0.227  0.001  0.041  0.132]\n"
     ]
    }
   ],
   "source": [
    "X = np.concatenate( (np.ones((X_train.shape[0],1)), X_train), axis=1)\n",
    "I = np.eye(p)\n",
    "I[0,0] = 0\n",
    "beta_R2 = np.linalg.solve(X.T@X + lam*I, X.T@Y_train)\n",
    "print(\"Estimation: \", beta_R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation:  [ 2.464  0.420  0.238 -0.047  0.162  0.227  0.001  0.041  0.132]\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cp\n",
    "\n",
    "lam = 24.25\n",
    "b = cp.Variable(p-1)\n",
    "b0 = cp.Variable(1)\n",
    "\n",
    "o = cp.Minimize(cp.sum_squares(b0+X_train@b-Y_train) + lam*cp.sum_squares(b))\n",
    "problem = cp.Problem(o)\n",
    "problem.solve(solver=cp.SCS,eps=1e-5)\n",
    "\n",
    "beta_R3 = np.concatenate((np.array(b0.value),b.value))\n",
    "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})\n",
    "print(\"Estimation: \",beta_R3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lasso regression\n",
    "It differs from the ridge regession by the use of a $\\ell_1$-norm regularization on the parameters\n",
    "$$\n",
    "\\min_{\\beta_0, \\, \\beta} \\sum_{i=1}^n \\left(y_i - \\beta_0 - x_i^\\top \\beta \\right)^2 + \\lambda \\sum_{j=1}^p |\\beta_j|\n",
    "$$\n",
    "An equivalent formulation (due to convexity of the fitting error term and the regularization) is:\n",
    "\\begin{eqnarray}\n",
    "\\min_{\\beta_0, \\, \\beta} & \\sum_{i=1}^n \\left(y_i - \\beta_0 - x_i^\\top \\beta \\right)^2 \\\\\n",
    "s.t. & \\sum_{j=1}^p |\\beta_j| \\leq t\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation:  [ 2.468  0.533  0.169 -0.000  0.002  0.094 -0.000 -0.000  0.000]\n"
     ]
    }
   ],
   "source": [
    "import cvxpy as cp\n",
    "\n",
    "Xn = (X_train - X_train.mean(axis=0))/X_train.std(axis=0)\n",
    "Yn = (Y_train - Y_train.mean())/Y_train.std()\n",
    "t = .7015\n",
    "n, p = Xn.shape\n",
    "\n",
    "b = cp.Variable(p)\n",
    "o = cp.Minimize(cp.sum_squares(Xn@b-Yn))\n",
    "c = [cp.norm(b, 1) <= t]\n",
    "problem = cp.Problem(o, c)\n",
    "problem.solve(solver=cp.SCS,eps=1e-5)\n",
    "\n",
    "ybar = Y_train.mean()\n",
    "stdy = Y_train.std()\n",
    "xbar = X_train.mean(0)\n",
    "stdx = X_train.std(0)\n",
    "b0_lasso =  ybar - stdy*((xbar/stdx).T@b.value)\n",
    "b_lasso = b.value*Y_train.std()/X_train.std(axis=0)\n",
    "beta_lasso = np.concatenate((np.array([b0_lasso]),b_lasso))\n",
    "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})\n",
    "print(\"Estimation: \",beta_lasso)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23.385463577394923\n",
      "Estimation:  [ 2.468  0.533  0.169  0.000  0.002  0.094  0.000 -0.000  0.000]\n"
     ]
    }
   ],
   "source": [
    "lam = c[0].dual_value\n",
    "print(lam)    # lam = 23.38\n",
    "o = cp.Minimize(cp.sum_squares(Xn@b-Yn) + lam*cp.norm(b, 1))\n",
    "problem = cp.Problem(o)\n",
    "problem.solve(solver=cp.SCS,eps=1e-5)\n",
    "\n",
    "b0_lasso =  ybar - stdy*((xbar/stdx).T@b.value)\n",
    "b_lasso = b.value*Y_train.std()/X_train.std(axis=0)\n",
    "beta_lasso = np.concatenate((np.array([b0_lasso]),b_lasso))\n",
    "np.set_printoptions(formatter={'float': '{: 0.3f}'.format})\n",
    "print(\"Estimation: \",beta_lasso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation:  [ 2.497  0.541  0.291 -0.153  0.214  0.318 -0.050  0.233 -0.061]\n"
     ]
    }
   ],
   "source": [
    "U, s, Vt = np.linalg.svd(Xn, full_matrices=False, compute_uv=True)\n",
    "Z = U*s\n",
    "theta = Z.T@Yn/np.diag(Z.T@Z)\n",
    "k = 7\n",
    "b_PCR = Vt.T[:,:k]@theta[:k]\n",
    "\n",
    "ybar = Y_train.mean()\n",
    "stdy = Y_train.std()\n",
    "xbar = X_train.mean(0)\n",
    "stdx = X_train.std(0)\n",
    "b0_PCR =  ybar - stdy*((xbar/stdx).T@b_PCR)\n",
    "b_PCR = b_PCR/X_train.std(axis=0)*stdy\n",
    "beta_PCR = np.concatenate((np.array([b0_PCR]),b_PCR))\n",
    "print(\"Estimation: \",beta_PCR)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PLS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.483  0.576  0.279 -0.179  0.201  0.299 -0.036  0.006  0.117]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cross_decomposition import PLSRegression\n",
    "\n",
    "pls2 = PLSRegression(n_components=3,scale=False)  # really better results with n_components=3\n",
    "#pls2.fit(Xn, Yn)\n",
    "X_train_1 = np.hstack((np.ones((n,1)),X_train))\n",
    "pls2.fit(X_train_1, Y_train)\n",
    "nt = X_test.shape[0]\n",
    "X_test_1 = np.hstack((np.ones((nt,1)),0*X_test))\n",
    "Y_pred = pls2.predict(X_test_1)\n",
    "\n",
    "pls2.coef_[0] = Y_pred[0]    # 2.452 - it doesn't fit the book !\n",
    "beta_pls = pls2.coef_.flatten()\n",
    "print(beta_pls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test_1 = np.hstack((np.ones((nt,1)),X_test))\n",
    "B = np.stack((b_ls.T , b_bs.T , beta_R.T , beta_lasso.T, beta_PCR.T, beta_pls.T)).T\n",
    "ytp = X_test_1@B\n",
    "err = ytp - np.outer(Y_test,np.ones(6))\n",
    "testErr = np.mean(err**2,axis=0)\n",
    "\n",
    "StdErr = np.std(err**2,axis=0, ddof=1)/np.sqrt(nt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "------------------------------------------------------------------------------------\n",
      "Term              LS      Best Subset    Ridge       Lasso        PCR         PLS       \n",
      "------------------------------------------------------------------------------------\n",
      "Intercept        2.465       2.477       2.452       2.468       2.497       2.483\n",
      "lcavol           0.680       0.740       0.420       0.533       0.541       0.576\n",
      "lweight          0.263       0.316       0.238       0.169       0.291       0.279\n",
      "age             -0.141       0.000      -0.047       0.000      -0.153      -0.179\n",
      "lbph             0.210       0.000       0.162       0.002       0.214       0.201\n",
      "svi              0.305       0.000       0.227       0.094       0.318       0.299\n",
      "lcp             -0.288       0.000       0.001       0.000      -0.050      -0.036\n",
      "gleason         -0.021       0.000       0.041      -0.000       0.233       0.006\n",
      "pgg45            0.267       0.000       0.132       0.000      -0.061       0.117\n",
      "------------------------------------------------------------------------------------\n",
      "Test Error       0.521       0.492       0.492       0.479       0.448       0.427\n",
      "Std Error        0.179       0.143       0.164       0.164       0.104       0.113\n",
      "------------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "dash = '-' * 84\n",
    "print(dash)\n",
    "print(\"{:<18s}{:<8s}{:<15s}{:<12s}{:<13s}{:<12s}{:<10s}\".format(\"Term\", \"LS\", \"Best Subset\", \"Ridge\", \"Lasso\", \"PCR\", \"PLS\"))\n",
    "print(dash)\n",
    "for k in range(variables.shape[0]+1):\n",
    "    if k==0:\n",
    "        print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(\"Intercept\", B[0,0],B[0,1],B[0,2],B[0,3],B[0,4],B[0,5]))\n",
    "    else:\n",
    "        print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(variables[k-1],B[k,0],B[k,1],B[k,2],B[k,3],B[k,4],B[k,5]))\n",
    "print(dash)\n",
    "print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(\"Test Error\",testErr[0],testErr[1],testErr[2],testErr[3],testErr[4],testErr[5]))\n",
    "print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(\"Std Error\",StdErr[0],StdErr[1],StdErr[2],StdErr[3],StdErr[4],StdErr[5]))\n",
    "print(dash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = \"https://i.stack.imgur.com/4Nooa.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improve the PLS results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.483  0.576  0.279 -0.179  0.201  0.299 -0.036  0.006  0.117]\n",
      "------------------------------------------------------------------------------------\n",
      "Term              LS      Best Subset    Ridge       Lasso        PCR         PLS       \n",
      "------------------------------------------------------------------------------------\n",
      "Intercept        2.465       2.477       2.452       2.468       2.497       2.483\n",
      "lcavol           0.680       0.740       0.420       0.533       0.541       0.576\n",
      "lweight          0.263       0.316       0.238       0.169       0.291       0.279\n",
      "age             -0.141       0.000      -0.047       0.000      -0.153      -0.179\n",
      "lbph             0.210       0.000       0.162       0.002       0.214       0.201\n",
      "svi              0.305       0.000       0.227       0.094       0.318       0.299\n",
      "lcp             -0.288       0.000       0.001       0.000      -0.050      -0.036\n",
      "gleason         -0.021       0.000       0.041      -0.000       0.233       0.006\n",
      "pgg45            0.267       0.000       0.132       0.000      -0.061       0.117\n",
      "------------------------------------------------------------------------------------\n",
      "Test Error       0.521       0.492       0.492       0.479       0.448       0.427\n",
      "Std Error        0.179       0.143       0.164       0.164       0.104       0.113\n",
      "------------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cross_decomposition import PLSRegression\n",
    "\n",
    "pls2 = PLSRegression(n_components=3,scale=False)  # really better results with n_components=3\n",
    "#pls2.fit(Xn, Yn)\n",
    "X_train_1 = np.hstack((np.ones((n,1)),X_train))\n",
    "pls2.fit(X_train_1, Y_train)\n",
    "nt = X_test.shape[0]\n",
    "X_test_1 = np.hstack((np.ones((nt,1)),0*X_test))\n",
    "Y_pred = pls2.predict(X_test_1)\n",
    "\n",
    "pls2.coef_[0] = Y_pred[0]    # 2.452\n",
    "beta_pls = pls2.coef_.flatten()\n",
    "print(beta_pls)\n",
    "\n",
    "X_test_1 = np.hstack((np.ones((nt,1)),X_test))\n",
    "B = np.stack((b_ls.T , b_bs.T , beta_R.T , beta_lasso.T, beta_PCR.T, beta_pls.T)).T\n",
    "ytp = X_test_1@B\n",
    "err = ytp - np.outer(Y_test,np.ones(6))\n",
    "testErr = np.mean(err**2,axis=0)\n",
    "\n",
    "StdErr = np.std(err**2,axis=0, ddof=1)/np.sqrt(nt)\n",
    "\n",
    "dash = '-' * 84\n",
    "print(dash)\n",
    "print(\"{:<18s}{:<8s}{:<15s}{:<12s}{:<13s}{:<12s}{:<10s}\".format(\"Term\", \"LS\", \"Best Subset\", \"Ridge\", \"Lasso\", \"PCR\", \"PLS\"))\n",
    "print(dash)\n",
    "for k in range(variables.shape[0]+1):\n",
    "    if k==0:\n",
    "        print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(\"Intercept\", B[0,0],B[0,1],B[0,2],B[0,3],B[0,4],B[0,5]))\n",
    "    else:\n",
    "        print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(variables[k-1],B[k,0],B[k,1],B[k,2],B[k,3],B[k,4],B[k,5]))\n",
    "print(dash)\n",
    "print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(\"Test Error\",testErr[0],testErr[1],testErr[2],testErr[3],testErr[4],testErr[5]))\n",
    "print(\"{:<10s}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}{:>12.3f}\".format(\"Std Error\",StdErr[0],StdErr[1],StdErr[2],StdErr[3],StdErr[4],StdErr[5]))\n",
    "print(dash)"
   ]
  },
  {
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}