practicals/week_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting started with the practicals\n",
    "\n",
    "***These notebooks are best viewed in Jupyter. GitHub might not display all content of the notebook properly.***\n",
    "\n",
    "## Goal of the practical exercises\n",
    "\n",
    "The exercises have two goals:\n",
    "\n",
    "1. Give you the opportunity to obtain 'hands-on' experience in implementing, training and evaluation machine learning models in Python. This experience will also help you better understand the theory covered during the lectures. \n",
    "\n",
    "2. Occasionally demonstrate some 'exam-style' questions that you can use as a reference when studying for the exam. Note however that the example questions are (as the name suggests) only examples and do not constitute a complete and sufficient list of 'things that you have to learn for the exam'. You can recognize example questions as (parts of) exercises by <font color=\"#770a0a\">this font color</font>.\n",
    "\n",
    "For each set of exercises (one Python notebook such as this one $==$ one set of exercises) you have to submit deliverables that will then be graded and constitute 25% of the final grade. Thus, the work that you do during the practicals has double contribution towards the final grade: as 30% direct contribution and as a preparation for the exam that will define the other 65% of the grade.\n",
    "\n",
    "## Deliverables\n",
    "\n",
    "For each set of exercises, you have to submit:\n",
    "1. Python functions and/or classes (`.py` files) that implement basic functionalities (e.g. a $k$-NN classifier) and \n",
    "2. A *single* Python notebook that contains the experiments, visualization and answer to the questions and math problems. *Do not submit your answers as Word or PDF documents (they will not be graded)*. The submitted code and notebook should run without errors and be able to fully reproduce the reported results.\n",
    "\n",
    "We recommend that you clone the provided notebooks (such as this one) and write your code in them. The following rubric will be used when grading the practical work:\n",
    "\n",
    "Component  | Insufficient | Satisfactory | Excellent\n",
    "--- | --- | --- | ---\n",
    "**Code** | Missing or incomplete code structure, runs with errors, lacks documentation | Self-contained, does not result in errors, contains some documentation, can be easily used to reproduce the reported results | User-friendly, well-structured (good separation of general functionality and experiments, i.e. between `.py` files and the Pyhthon notebook), detailed documentation, optimized for speed, use of a version control system (such as GitHub)\n",
    "**Answers to questions** | Incorrect, does not convey understanding of the material, appears to be copied from another source | Correct, conveys good understanding of the material, description in own words | Correct, conveys excellent level of understanding, makes connections between topics\n",
    "\n",
    "## A word on notation\n",
    "\n",
    "When we refer to Python variables, we will use a monospace font. For example, `X` is a Python variable that contains the data matrix. When we refer to mathematical variables, we will use the de-facto standard notation: $a$ or $\\lambda$ is a scalar variable, $\\boldsymbol{\\mathrm{w}}$ is a vector and $\\boldsymbol{\\mathrm{X}}$ is a matrix (e.g. a data matrix from the example above). You should use the same notation when writing your answers and solutions.\n",
    "\n",
    "# Two simple machine learning models\n",
    "\n",
    "## Preliminaries\n",
    "\n",
    "Throughout the practical curriculum of this course, we will use the Python programming language and its ecosystem of libraries for scientific computing (such as `numpy`, `scipy`, `matplotlib`, `scikit-learn` etc). The practicals for the deep learning part of the course will use the `keras` deep learning framework. If you are not sufficiently familiar with this programming language and/or the listed libraries and packages, you are strongly advised to go over the corresponding tutorials from the ['Essential skills'](https://github.com/tueimage/essential-skills) module (the `scikit-learn` library is not covered by the tutorial, however, an extensive documentation is available [here](https://scikit-learn.org/stable/documentation.html).\n",
    "\n",
    "In this first set of exercises, we will use two toy datasets that ship together with `scikit-learn`. \n",
    "\n",
    "The first dataset is named `diabetes` and contains 442 patients described with 10 features: age, sex, body mass index, average blood pressure, and six blood serum measurements. The target variable is a continuous quantitative measure of the disease (diabetes) progression one year after the baseline measurements were recorded. More information is available [here](https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/datasets/descr/diabetes.rst) and [here](https://www4.stat.ncsu.edu/~boos/var.select/diabetes.html).\n",
    "\n",
    "The second dataset is named `breast_cancer` and is a copy of the UCI ML Breast Cancer Wisconsin (Diagnostic) datasets (more infortmation is available [here](https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/datasets/descr/breast_cancer.rst) and [here](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic)). The datasets contains of 569 instances represented with 30 features that are computed from a images of a fine needle aspirate of a breast mass. The features describe characteristics of the cell nuclei present in the image. Each instance is associated with a binary target variable ('malignant' or 'benign'). \n",
    "\n",
    "You can load the two datasets in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import load_diabetes, load_breast_cancer\n",
    "\n",
    "diabetes = load_diabetes()\n",
    "\n",
    "breast_cancer = load_breast_cancer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the majority of the exercises in this course, we will use higher-level libraries and packages such as `scikit-learn` and `keras` to implement, train and evaluate machine learning models. However, the goal of this first set of exercises is to illustrate basic mathematical tools and machine learning concepts. Because of this, we will impose a restriction of only using basic `numpy` functionality. Furthermore, you should as much as possible restrict the use of for-loops (e.g. use a vector-to-matrix product instead of a for loop when appropriate).\n",
    "\n",
    "If `X` is a 2D data matrix, we will use the convention that the rows of the matrix contain the samples (or instances) and the columns contain the features (inputs to the model). That means that a data matrix with a shape `(122, 13)` represents a dataset with 122 samples, each represented with 13 features. Similarly, if `Y` is a 2D matrix containing the targets, the rows correspond to the samples and the columns to the different targets (outputs of the model). Thus, if the shape of `Y` is `(122, 3)` that means that there are 122 samples and each sample is has 3 targets (note that in the majority of the examples we will only have a single target and thus the number of columns of `Y` will be 1).\n",
    "\n",
    "You can obtain the data and target matrices from the two datasets in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(442, 10)\n",
      "(442, 1)\n"
     ]
    }
   ],
   "source": [
    "X = diabetes.data\n",
    "Y = diabetes.target[:, np.newaxis]\n",
    "\n",
    "print(X.shape)\n",
    "print(Y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to only use a subset of the available features, you can obtain a reduced data matrix in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(442, 1)\n",
      "(442, 2)\n"
     ]
    }
   ],
   "source": [
    "# use only the fourth feature\n",
    "X = diabetes.data[:, np.newaxis, 3]\n",
    "print(X.shape)\n",
    "\n",
    "# use the third, and tenth features\n",
    "X = diabetes.data[:, (3,9)]\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Question***: Why we need to use the `np.newaxis` expression in the examples above? \n",
    "\n",
    "Answer: `np.newaxis` is used to ensure that X is a 2D matrix instead of a 1D vector, and that this dimension has the correct size, which in this case is the default 1.\n",
    "\n",
    "Note that in all your experiments in the exercises, you should use and independent training and testing sets. You can split the dataset into a training and testing subsets in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(300, 1)\n",
      "(300, 1)\n",
      "(142, 1)\n",
      "(142, 1)\n"
     ]
    }
   ],
   "source": [
    "# use the fourth feature\n",
    "# use the first 300 training samples for training, and the rest for testing\n",
    "X_train = diabetes.data[:300, np.newaxis, 3]\n",
    "y_train = diabetes.target[:300, np.newaxis]\n",
    "X_test = diabetes.data[300:, np.newaxis, 3]\n",
    "y_test = diabetes.target[300:, np.newaxis]\n",
    "print(X_train.shape)\n",
    "print(y_train.shape)\n",
    "print(X_test.shape)\n",
    "print(y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "### Linear regression\n",
    "\n",
    "Implement training and evaluation of a linear regression model on the diabetes dataset using only matrix multiplication, inversion and transpose operations. Report the mean squared error of the model.\n",
    "\n",
    "To get you started we have implemented the first part of this exercise (fitting of the model) as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 152.34786452]\n",
      " [ -16.57607993]\n",
      " [-254.66532396]\n",
      " [ 560.98630022]\n",
      " [ 278.91811152]\n",
      " [-393.41357305]\n",
      " [  97.05460405]\n",
      " [ -19.0023093 ]\n",
      " [ 169.46450327]\n",
      " [ 632.95050374]\n",
      " [ 114.21638941]]\n",
      "Mean Square Error: 2794.5690145007343\n"
     ]
    }
   ],
   "source": [
    "# add subfolder that contains all the function implementations\n",
    "# to the system path so we can import them\n",
    "import sys\n",
    "sys.path.append('code/')\n",
    "\n",
    "# the actual implementation is in linear_regression.py,\n",
    "# here we will just use it to fit a model\n",
    "from linear_regression import *\n",
    "\n",
    "# load the dataset\n",
    "# same as before, but now we use all features\n",
    "X_train = diabetes.data[:300, :]\n",
    "y_train = diabetes.target[:300, np.newaxis]\n",
    "X_test = diabetes.data[300:, :]\n",
    "y_lest = diabetes.target[300:, np.newaxis]\n",
    "\n",
    "w = lsq(X_train, y_train)\n",
    "\n",
    "# print the parameters\n",
    "print(w)\n",
    "\n",
    "y_predicted = predict(X_test,w)\n",
    "meanSquareError = msq(y_test,y_predicted)\n",
    "print(\"Mean Square Error: \"+str(meanSquareError))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weighted linear regression\n",
    "\n",
    "Assume that in the dataset that you use to train a linear regression model, there are identical versions of some samples. This problem can be reformulated to a weighted linear regression problem where the matrices $\\boldsymbol{\\mathrm{X}}$ and $\\boldsymbol{\\mathrm{Y}}$ (or the vector $\\boldsymbol{\\mathrm{y}}$ if there is only a single target/output variable) contain only the unique data samples, and a vector $\\boldsymbol{\\mathrm{d}}$ is introduced that gives more weight to samples that appear multiple times in the original dataset (for example, the sample that appears 3 times has a corresponding weight of 3). \n",
    "\n",
    "<p><font color='#770a0a'>Derive the expression for the least-squares solution of a weighted linear regression model (note that in addition to the matrices $\\boldsymbol{\\mathrm{X}}$ and $\\boldsymbol{\\mathrm{Y}}$, the solution should include a vector of weights $\\boldsymbol{\\mathrm{d}}$).</font></p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to solve the weighted mean square error formula:\n",
    "$$WMSE(\\beta,d1,...,dn) = \\frac{1}{n}\\sum_{i=1}^nd_i(y_i-\\boldsymbol{\\mathrm{x_i}}\\boldsymbol{\\mathrm{\\beta}})$$\n",
    "Matrix $\\boldsymbol{\\mathrm{W}}$ will be the matrix with the values of vector $\\boldsymbol{\\mathrm{d}}$ on the diagonal. This gives:\n",
    "$$WMSE = n^{-1}(\\boldsymbol{\\mathrm{y}}-\\boldsymbol{\\mathrm{X}}\\boldsymbol{\\mathrm{\\beta}})^T\\boldsymbol{\\mathrm{W}}(\\boldsymbol{\\mathrm{y}}-\\boldsymbol{\\mathrm{X}}\\boldsymbol{\\mathrm{\\beta}})$$\n",
    "$$=\\frac{1}{n}(\\boldsymbol{\\mathrm{y}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{y}}-\\boldsymbol{\\mathrm{y}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{X}}\\boldsymbol{\\mathrm{\\beta}}-\\boldsymbol{\\mathrm{\\beta}}^T\\boldsymbol{\\mathrm{X}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{y}}+\\boldsymbol{\\mathrm{\\beta}}^T\\boldsymbol{\\mathrm{X}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{X}}\\boldsymbol{\\mathrm{\\beta}}$$\n",
    "Differentiate with respect to $\\boldsymbol{\\mathrm{\\beta}}$ to get the gradient:\n",
    "$$\\nabla_\\boldsymbol{\\mathrm{\\beta}}WMSE = \\frac{2}{n}(-\\boldsymbol{\\mathrm{X}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{y}}+\\boldsymbol{\\mathrm{X}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{X}}\\boldsymbol{\\mathrm{\\beta}})$$\n",
    "Then it follows that for the optimum at zero $\\boldsymbol{\\mathrm{\\beta}}$ is:\n",
    "$$\\boldsymbol{\\mathrm{\\hat{\\beta}}} = (\\boldsymbol{\\mathrm{X}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{X}})^{-1}\\boldsymbol{\\mathrm{X}}^T\\boldsymbol{\\mathrm{W}}\\boldsymbol{\\mathrm{y}}$$\n",
    "Or when $\\boldsymbol{\\mathrm{d}}$ must truly appear:\n",
    "$$\\boldsymbol{\\mathrm{\\hat{\\beta}}} = (\\boldsymbol{\\mathrm{X}}^Tdiag(\\boldsymbol{\\mathrm{d}})\\boldsymbol{\\mathrm{X}})^{-1}\\boldsymbol{\\mathrm{X}}^Tdiag(\\boldsymbol{\\mathrm{d}})\\boldsymbol{\\mathrm{y}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $k$-NN classification\n",
    "\n",
    "Implement a $k$-Nearest neighbors classifier from scratch in Python using only basic matrix operations with `numpy` and `scipy`. Train and evaluate the classifier on the breast cancer dataset, using all features. Show the performance of the classifier for different values of $k$ (plot the results in a graph). Note that for optimal results, you should normalize the features (e.g. to the $[0, 1]$ range or to have a zero mean and unit standard deviation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8lfWZ9/HPNwlJICQsWYCwrwkBFStqXUlsp9Wqdev00XG6L0+ndZ7OjHaqdapTOo7tTJeZLtNn7FRbOnVsH6czWodWKU3ErRYoqAQSQBaBAAlbQoCQ7Xr+OHfwmAZyQnJyn+V6v17nlXPv1w3Juc5vuX8/mRnOOefc2coIOwDnnHPJzROJc865QfFE4pxzblA8kTjnnBsUTyTOOecGxROJc865QfFE4lwCUcQjkg5L+l3Y8TgXC08kLiVIqgk+fHPCjmWQLgf+CJhiZhcN9mSSZkgySVnBsiR9W1KdpMmSPhxs/1yv43ZLqgze/22wzx9Hbc8K1s0YbIwu+XkicUkv+DC7AjDgvcN87awhPuV0YIeZHRvqWCQJ+FegElhiZnuCTYeAz0sqOMPhh4ClkjIHGpdLfZ5IXCr4IPBb4IfAh6I3SBop6euSdkpqlvS8pJHBtsslvSjpiKRdkj4crK+R9PGoc3xY0vNRyybpM5K2AFuCdf8cnKNF0lpJVwTrJ0o6Lqkw6vgLJDVJGtEr1o8B/wZcIqlV0peC9Z+QtFXSIUlPSio9UyynkRn8+ywGKs1sf9S2TcBLwF+e4fhfAe3An55hH5emPJG4VPBB4CfB692SJkRt+xpwAXApMB74a6Bb0jTgl8C3gWJgEbB+ANe8EbgYqAiWVwfnGA88Cvw/Sblmtg+oAd4fdeyfAo+ZWUf0Cc3sB8CngJfMbLSZ3S/pKuDB4PhJwE7gsX5i6ctPgHLgKjM72Mf2LwJ/KWn8aY63YJ/7eydA5zyRuKQm6XIi1UE/M7O1wOvAnwTbMoCPAp81sz1m1mVmL5rZSeB24Ndm9h9m1mFmB81sIInkQTM7ZGYnAMzs34NzdJrZ14EcoCzY90cE3+SDqqHbgB/HeJ3bgYfN7PdB3PcQKbHMOF0sp/EuIv9GR/raGNz7M8DnT3cCM3sSaAI+frp9XHryROKS3YeAZ8zsQLD8KG9WbxUBuUSSS29TT7M+VruiFyTdKWlTUH12BBgTXB/gCaBC0iwiDenNZhZrj6xSIqUQAMysFTgITD5dLKdxHZHSxEfPsM99wJ9JmniGff4GuJfIv6tzAAx1Q6FzwyZo63g/kClpX7A6Bxgr6TzgNaANmA280uvwXcDpekUdA0ZFLff1wXpq2OygPeTzwDuAWjPrlnQYEICZtUn6GZHSRTmxl0YAGoiUuHqulQcUAnui9ollCO8XgeuBpyS1mdmjf3BDZnWSfg584XQnMbMVkrYCn44xfpcGvETiktmNQBeRtoFFwWs+8BzwQTPrBh4GviGpVFKmpEuCLsI/Ad4p6f1BV9ZCSYuC864HbpY0StIc4GP9xJEPdBKp9smSdB/QuwfUMuDDRHqV/fsA7vFR4COSFgVx/z3wspntGMA5ADCzZ4GbgYckve80u30J+Agw9gynupdIW5NzgCcSl9w+BDxiZm+Y2b6eF/Ad4PagO+xdREomq4l0Yf0qkGFmbwDvAe4M1q8HzgvO+00iPZT2E2nf+Ek/cTxNpOF+M5FqqDZ6VTeZ2QtAN/D7gSQBM1tJpJH7P4G9REpXt8Z6fB/nWwH8L+CHkq7vY/t2IiWmvDOc4wXAH5Z0p8gntnJueEj6DfComf1b2LE4N5Q8kTg3DCRdCKwApprZ0bDjcW4oedWWc3Em6UfAr4G/8CTiUpGXSJxzzg2Kl0icc84NSlo8R1JUVGQzZswIOwznnEsqa9euPWBmxf3tlxaJZMaMGaxZsybsMJxzLqlI2tn/Xl615ZxzbpA8kTjnnBsUTyTOOecGxROJc865QfFE4pxzblA8kTjnnBsUTyTOOecGJS2eI3HwP6/upX5fy4COGZ2bxUcvm0lWZnJ/39jW1Mrm/a1cvfBME/85586WJ5I08NLrB/nMo78HQIrtmJ4h2GYU5vGuBcn9AfzVX9XxzMb9PPGZyzh3ypnma3LOnQ1PJCmuq9v40i9qmTx2JCvvXELuiMyYjmvv7Ob8pc9Qs7kpqRNJe2c3z285gBl86RcbefxTl6BYs6lzLibJXWfh+vXY6jeo23eUL7xnfsxJBCA7K4PL5xZRU9dIMo8QvWbHIY61d/Guigms3XmYJ19pCDsk51KOJ5IU1nyig68/s5mLZo7nPecMvFRRVVZCQ3Mbm/e3xiG64VFd30h2ZgZff/95LJxcwFd+Wcfx9s6ww3IupXgiSWHfWrmFw8fbuf/6irOqzqksKwEiH8bJqqa+iYtnjSc/dwT3X7+Avc1t/Ouz28IOy7mU4okkRW1tbOVHL+7g1gunsaB0zFmdY+KYXOZPKqAmSRPJrkPH2dLYypJ5kVGwL5wxnuvOncT/ffZ19hw5EXJ0zqWOuCYSSVdLqpe0VdLdfWyfLmmlpFcl1UiaErVtmqRnJG2StFHSjGD9DyVtl7Q+eC2K5z0kq7/7n42MHJHJne+aN6jzVJUVs2bHYVraOoYosuFTs7kJgKryklPr7nnPfAAeXL4plJicS0VxSySSMoHvAtcAFcBtkip67fY1YJmZnQssBR6M2rYM+Eczmw9cBER/Lf6cmS0KXuvjdQ/JqrqukZr6Jj77zrkUjc4Z1Lkqy0ro7DZe2HJgiKIbPjV1jUwbP4pZRXmn1k0eO5JPLZnNU6/u5XfbD4UYnXOpI54lkouArWa2zczagceAG3rtUwGsDN5X92wPEk6Wma0AMLNWMzsex1hTRntnN19+aiOzivL44CUzBn2+t00bS35uVtK1k7R1dPHC6weoKiv+g/ahTy2ZzaQxuXzpF7V0dSdvjzTnEkU8E8lkYFfU8u5gXbRXgFuC9zcB+ZIKgXnAEUk/l7RO0j8GJZweDwTVYd+U1OdXbkmflLRG0pqmpqahuaMksOylHWw7cIwvXldBdtbg/3uzMjO4cl4x1fVNSdUN+OXth2jr6KYyqlqrx8jsTO55z3xqG1p4fO2uPo52zg1EPBNJX92Een8S3QUskbQOWALsATqJPCh5RbD9QmAW8OHgmHuA8mD9eODzfV3czB4ys8Vmtri4uN8ph1PCgdaT/PPKLSyZV/yWdoHBqioroenoSWobBjbESpiq6xrJycrgklmFfW6//txJLJ4+jn98uj4p23+cSyTxTCS7galRy1OAtzwNZmYNZnazmZ0P3Busaw6OXRdUi3UC/w28Ldi+1yJOAo8QqUJzwNef2cyJ9i6+eN38IT1vT6+nZOq9VVPfyKWzC0/7EKYk7r9+AQePtfOd32wd5uicSy3xTCSrgbmSZkrKBm4FnozeQVKRpJ4Y7gEejjp2nKSeosRVwMbgmEnBTwE3AhvieA9Jo7ahmcdWv8EHL5nBnJL8IT13cX4O504ZQ3V9clQRbj9wjB0Hj/dbKjtnyhj++IIpPPLCdrY1Je9Dl86FLW6JJChJ3AE8DWwCfmZmtZKWSnpvsFslUC9pMzABeCA4totItdZKSa8RqSb7fnDMT4J1rwFFwN/F6x6ShZnxpV9sZNyobD77jrlxuUZlWQnr3jjM4WPtcTn/UKqui5ScKuf1X71317vLyMnK5IH/8e7Azp2tuA7aaGbLgeW91t0X9f5x4PHTHLsCOLeP9VcNcZhJb/lr+/jd9kM8cNNCxowaEZdrVJUV862VW1i1pYkbFvXuM5FYajY3Mas4j2mFo/rdtyQ/lz+/ag4P/rKOZzc3narGc87Fzp9sT3JtHV38/fJNlE/M59YLp8XtOudOGcv4vGxqErx663h7J7/ddpCqstg7G3z4shnMKBzFl5/aSEdXdxyjcy41eSJJct9ftY09R05w3/UVZGbEb3j0zAyxZF4xz25uojuBn7146fWDtHd2DyiR5GRlcu+1FWxtbOXff7szjtE5l5o8kSSxvc0n+Jea17lm4UQunV0U9+tVlhVz6Fg7r+5pjvu1zlZ1fSOjsjO5cOa4AR33zvklXDG3iG+u2MyhJGgHci6ReCJJYl/9ZR1dZnzhPUPb3fd0rpxbTIbebMxONGZGdV0Tl80pIicr9rlXINId+IvXVXCsvYtvrKiPU4TOpSZPJElq7c5D/Pf6Bj5xxUymju+/UXkojMvLZtHUsQn7PMnWxlb2HDkxoGqtaPMm5POBt0/n0ZffoG6A89s7l848kSSh7u5Id9+S/Bw+XTlnWK9dVVbCK7ubaTp6clivG4ue8cAqy86+59VfvHMuBSNH8KUnNybVkDDOhcnnbE9CP1+3h1d3N/ON959HXs7w/hdWlZfw9RWbWbW5iVsumNL/AcOouq6J8on5lI4dedbnGDsqm7/6o3nc90Qtf/34q4zLyx7CCIfGe88rZeHks5tjxrl48ESShH7y8k7KJ+ZzYwjPc1RMKqA4P4fq+saESiRH2zpYveMQH79i1qDP9ScXTePXmxp56tW9QxDZ0Oro6uaJ9Xv4zZ2Vw/4lwrnT8d/EJHPoWDvrdx3hs++YS0Ycu/ueTkaGqJxXzNO1++js6iYrMzFqR1/YepDObqNqENVaPbIyM1j20cQcwm3tzsPc8r0X+V7N69z17rKww3EO8DaSpLNqcxNmnHWD8lCoKi+hpa2TdbuOhBZDbzX1jeTnZPG26QPr9ptsLpg+jhsXlfLQc9vYdcin6HGJwRNJkqmub6QwL5tzQqwjv3xuEZkZSphuwGZGdX0jV8wrYkSClJDi6fPXlJMp8fc+XbBLEKn/V5dCurotMh5UWXEo1Vo9CnJHsHj6uIQZDXjT3qPsbzlJZYiltOE0acxIPl05m19u2MdLrx8MOxznPJEkk/W7jnDkeEeo1Vo9qspL2LS3hX3NbWGH8ma33zQacPETV85i8tiRPl2wSwieSJJITX0jGYo8YR62nmSWCA8n1tQ3snByASUFuWGHMmxyR2Ry77Xzqdt3lMdWvxF2OC7NeSJJItX1jVwwfVzchoofiHkTRlM6Jjf00YCbj3ewdufhhCilDbdrFk7k4pnj+drT9TQf9+mCXXg8kSSJxpY2NuxpSZh2AElUlpfw/NYDtHeGN/T6qi1NdBsJ8+8ynCRx3/UVHDnRwT+v3BJ2OC6NeSJJEjWbI9/8BzP8x1CrnFdM68lO1uw8FFoM1fWNjB01gkVTx4YWQ5gWlI7h1gunseylHWxtPBp2OC5NeSJJEs/WN1GSn0PFpIKwQznlsjlFjMhUaNVb3d3GqmBWw3jOxZLo7nrXPEZmZ/Llp7w7sAuHJ5Ik0NHVzaotTVSVlSAlzgdmXk4WF88sDO15kg0NzRxobU/L9pFohaNz+Ow75vLs5qaEebbHpRdPJEng9zsPc7Stk6ryxKnW6lFZVsyWxtZQnrKurmtCgivTqNvv6XzwkhnMKs7jy09tDLXNyqUnTyRJoLq+iawMcdmc+M+COFBV5UE34M3DX71VXd/IecFc8ukuOyuDL15bwbYDx1j20o6ww3FpJq6JRNLVkuolbZV0dx/bp0taKelVSTWSpkRtmybpGUmbJG2UNCNYP1PSy5K2SPqppJT/FKmpb+TCGePJzw2/229vs4rymDZ+FDXDXKVysPUkr+w+kvbVWtGqykuoLCvmn3+9hQOtiTdfjEtdcUskkjKB7wLXABXAbZIqeu32NWCZmZ0LLAUejNq2DPhHM5sPXAT0fFJ9Ffimmc0FDgMfi9c9JIKGIyeo23c0Iau1INIFtaqsmBdeP0BbR9ewXXfVlmDwygT9dwnL31xbwYmOLr7+jE8X7IZPPEskFwFbzWybmbUDjwE39NqnAlgZvK/u2R4knCwzWwFgZq1mdlyRluargMeDY34E3BjHewhdT4+oRP7mXVleQltHNy9vH75uwNV1TRSNzmZhqU/wFG1OyWg+dOkMHlu9iw17msMOx6WJeCaSycCuqOXdwbporwC3BO9vAvIlFQLzgCOSfi5pnaR/DEo4hcARM+s8wzkBkPRJSWskrWlqSozBBc9GdX0jk8eOZE7J6LBDOa1LZhWSk5UxbD2GTg1eOa8k1MErE9X/ecdcxo3KZulTPl2wGx7xTCR9/YX3/q2+C1giaR2wBNgDdBKZcOuKYPuFwCzgwzGeM7LS7CEzW2xmi4uLk7P642RnFy9sPUBVeXFCdfvtLXdEJpfOLhy2cbfW7zpM84kOr9Y6jTEjR3Dnu+bxu+2HWP7avrDDcWkgnolkNzA1ankK0BC9g5k1mNnNZnY+cG+wrjk4dl1QLdYJ/DfwNuAAMFZS1unOmUpWbz/M8fauhK7W6lFVXsKOg8fZfuBY3K9VXddEZoa4Yo4nktO59cJplE/M5++XbxrWtiuXnuKZSFYDc4NeVtnArcCT0TtIKpLUE8M9wMNRx46T1PNJcRWw0SLl9GrgfcH6DwFPxPEeQlVT30h2VgaXzC4MO5R+Vc4bvtGAazY3csG0xBi8MlFlZoj7r1/AniMneGjVtrDDcSkubnO2m1mnpDuAp4FM4GEzq5W0FFhjZk8ClcCDkgxYBXwmOLZL0l3AyqCBfS3w/eDUnwcek/R3wDrgB/G6h7BV1zdy8czxjMqO23/TkJlWOIpZxXn8+KWdcS2VdJuxYU8Ln/P5yvt1yexC3nPORP6lZivvXzyViWPSZ5h9N7zi+gllZsuB5b3W3Rf1/nHe7IHV+9gVwLl9rN9GpEdYSnvj4HFebzrG7RdPDzuUmN1+8XS+85st/OKV+NY2Th47kmvPmRTXa6SKO99VxvLX9rFi4z4+cMmMsMNxKSrxv+qmqZrNkSqinifHk8HHLp/Jxy6fGXYYLsqsojzGjRpBbUNL2KG4FOZDpCSo6rpGZhSOYmZRXtihuCQmiQWlY9jQ4M+UuPjxRJKA2jq6ePH1g2k5WZMbegtKC9i8r5WOLh/M0cWHJ5IE9NK2g5zs7E6qai2XuCpKC2jv6mbL/tawQ3EpyhNJAqqpayR3RAYXzxwfdiguBSycHBlGxqu3XLx4IkkwZkZ1fROXzS4id0Rm2OG4FDCzMI9R2Zls9AZ3FyeeSBLMtgPHeOPQcSq9WssNkYwMMX9SAbVeInFx4okkwfSM9lvps/65IbSgtICNDS10d/sgjm7oeSJJMDX1jcwtGc3U8aPCDsWlkIWlYzjW3sWOg/EfC82lH08kCeTYyU5e3nbIe2u5IVdRWgDgDya6uPBEkkBefP0g7V3dXq3lhty8CfmMyJQnEhcXnkgSSHV9I3nZmSye4d1+3dDKzspg3oR8b3B3ceGJJEGYGTV1jVw+t4jsLP9vcUNvQWkBtQ0tPmuiG3L+iZUgNu9vpaG5LSkmsXLJaUHpGA4da2dfS1vYobgU44kkQVQHE0L5+FouXhb0NLjv8XYSN7Q8kSSI6rpG5k8q8MmHXNzMn1SA5EOluKHniSQBtLR1sHbnYarKvLeWi5+8nCxmFuV5zy035DyRJIAXthygs9v8+REXdwtKx/iYW27IeSJJANX1jRTkZnH+1LFhh+JS3MLSAvYcOcHhY+1hh+JSiCeSkPWM9nvlvGKyMv2/w8XXgtLIkPJeveWGUlw/uSRdLale0lZJd/exfbqklZJelVQjaUrUti5J64PXk1Hrfyhpe9S2RfG8h3irbWih6ehJ7/brhsWpnlve4O6GUFa8TiwpE/gu8EfAbmC1pCfNbGPUbl8DlpnZjyRdBTwIfCDYdsLMTpckPmdmj8cr9uG0eschAC6bUxRyJC4djMvLpnRMrpdI3JCKZ4nkImCrmW0zs3bgMeCGXvtUACuD99V9bE95tQ0tFI3O8W6/btgsmDzGuwC7IRXPRDIZ2BW1vDtYF+0V4Jbg/U1AvqTCYDlX0hpJv5V0Y6/jHgiqw74pKaevi0v6ZHD8mqampkHeSvzUNrScqm5wbjgsKC1g+4FjHDvZGXYoLkXEM5Goj3W9B/m5C1giaR2wBNgD9Px2TzOzxcCfAP8kaXaw/h6gHLgQGA98vq+Lm9lDZrbYzBYXFyfm8xknO7vYsv+oJxI3rBaUjsEM6vZ59ZYbGvFMJLuBqVHLU4CG6B3MrMHMbjaz84F7g3XNPduCn9uAGuD8YHmvRZwEHiFShZaUNu9rpbPbWDh5TNihuDSycLLPTeKGVjwTyWpgrqSZkrKBW4Eno3eQVCSpJ4Z7gIeD9eN6qqwkFQGXARuD5UnBTwE3AhvieA9x1dNzxkskbjhNLMhlfF42G/Z4O4kbGnHrtWVmnZLuAJ4GMoGHzaxW0lJgjZk9CVQCD0oyYBXwmeDw+cC/Suomkuy+EtXb6yeSiolUna0HPhWve4i3DQ3N5OdkMXWcT6vrho+kU0PKOzcU4pZIAMxsObC817r7ot4/DvxBN14zexE45zTnvGqIwwxNbUML80sLyMjoqznJufipKC3g4ee3097Z7fPfuEHz36CQdHUbdXuPsrDU20fc8FtYOoaOLmNL49GwQ3EpwBNJSLYfaOVER5e3j7hQ+Nwkbih5IgnJhuAPeMFkTyRu+M0ozCMvO9OHSnFDwhNJSGobmsnJymBO8eiwQ3FpKCNDzJ/kDe5uaHgiCUltQwvlE/N9xF8XmoWTx7Bxbwtd3b2fE3ZuYGL6FJP0n5KujXrmww2CmVHb0EKFN7S7EFWUFnC8vYsdB4+FHYpLcrEmhu8RGapki6SvSCqPY0wpb/fhEzSf6PCGdheqN4eU9+otNzgxJRIz+7WZ3Q68DdgBrJD0oqSPSBoRzwBTUc8frg+N4sI0tySf7MwMb3B3gxZzVVUwKu+HgY8D64B/JpJYVsQlshS2saGZzAxRPjE/7FBcGsvOymDexNHeBdgNWkxPtkv6OZERd38MXG9me4NNP5W0Jl7BpaoNDS3MLs4jd0Rm2KG4NLdg0hie2bgPMyMyfJ1zAxdrieQ7ZlZhZg9GJREAgqHe3QDUNjSfmjvbuTAtmFzA4eMd7G1uCzsUl8RiTSTzJY3tWQhG5/10nGJKaQdaT7K/5aQ3tLuE0POFxhvc3WDEmkg+YWZHehbM7DDwifiElNp6/mC9ROISwfxJ+Uj4kPJuUGJNJBmKqkCVlAlkxyek1NbzB1vhJRKXAEZlZzGrKC8pSiTbDxzjGys20+0PUCacWBPJ08DPJL1D0lXAfwC/il9YqWtjQwvTxo9izEjvNe0Sw8LJY9iYBF2Af/D8Nr61cguveukp4cSaSD4P/Ab4MyKTT60E/jpeQaWySEO7l0Zc4lhQWkBDcxuHjrWHHcppmRk19U0AVNc1hhyN6y3WBxK7zex7ZvY+M7vFzP7VzLriHVyqaWnrYMfB455IXEJ5s8E9cb/pv97Uyu7DJ5CgZnNT2OG4XmIda2uupMclbZS0recV7+BSzSZvaHcJKBmGSqmuiySPWy+cyqu7j3Cg9WTIEblosVZtPUJkvK1OoApYRuThRDcAp3ps+RwkLoGMHZXN5LEjEzuR1DdSNiGf2y+ejhms8lJJQok1kYw0s5WAzGynmf0tkDJzpw+X2oYWivNzKMnPDTsU595iQWkBtQnaiH20rYPVOw5RWV5MxaQCivNzqK73RJJIYk0kbcEQ8lsk3SHpJqAkjnGlJG9od4lqQekYth88xrGTnWGH8gde2HqQji6jqqyEjAxROa+YVZub6OzqDjs0F4g1kfwFMAr4P8AFwJ8CH+rvIElXS6qXtFXS3X1sny5ppaRXJdVImhK1rUvS+uD1ZNT6mZJelrRF0k8lJcXzLG0dXWxpbPVE4hLSwskFmMGmvYlXvVVT30h+ThYXTB8HQFV5Cc0nOli/60g/R7rh0m8iCR4+fL+ZtZrZbjP7SNBz67cxHPdd4BqgArhNUkWv3b4GLDOzc4GlwINR206Y2aLg9d6o9V8Fvmlmc4HDwMf6u4dEsHn/Ubq6jYXe0O4SUE8HkER7wt3MqK5v5Ip5RYwIZhO9fG4RmRmiut67ASeKfhNJ0M33gugn22N0EbDVzLaZWTvwGHBDr30qiDyTAlDdx/a3CGK4Cng8WPUj4MYBxhUKHxrFJbIJBTkU5mUnXIP7pr1H2d9yksqyN2vSC3JHsHj6uFM9uVz4Yq3aWgc8IekDkm7uefVzzGRgV9Ty7mBdtFeAW4L3NwH5wbwnALmS1kj6raSeZFEIHDGznorcvs4JgKRPBsevaWoK/xduw55m8nOzmDp+ZNihOPcHJFFRWpBwiaSn1FE5r/gt6yvLSti4t4V9PmpxQog1kYwHDhIpDVwfvK7r55i+SjC9B8m5C1giaR2wBNhDpIsxwLRgiPo/Af5J0uwYzxlZafaQmS02s8XFxcV97TKsahtaWFBa4HM+uIS1cPIYtjQe5WRn4jxr/Gx9EwtKCygpeGtPx6ryyN/0s5u9eisRxDSxlZl95CzOvRuYGrU8BWjodd4G4GYASaOBW8ysOWobZrZNUg1wPvCfwFhJWUGp5A/OmYg6u7qp29fC7RdPDzsU505rQWkBHV3Glv2tCTENdPPxDta+cZg/WzL7D7aVTchn0phcquua+F8XTgshOhct1hkSH6GPb/5m9tEzHLYamCtpJpGSxq1EShfR5y0CDplZN3AP8HCwfhxw3MxOBvtcBvyDmZmkauB9RNpcPgQ8Ecs9hGnbgWO0dXR7jy2X0KKHSkmERPLc1ia6uu1U6SOaJCrLSvjFKw20d3aTnRXzrOEuDmL9138K+J/gtRIoAFrPdEBQYriDyMjBm4CfmVmtpKWSenphVQL1kjYDE4AHgvXzgTWSXiHSCP8VM9sYbPs88FeSthJpM/lBjPcQmp4xjLyh3SWy6eNHMTonK2HaSarrmhg7agSLpo7rc3tVWTGtJztZu/PwMEfmeou1aus/o5cl/Qfw6xiOWw4s77Xuvqj3j/NmD6zofV4EzjnNObcR6RGWNGr3tJCTlcHs4rywQ3HutDIyRMWkxGhw7+42nt3cyJVzi8nM6Ltd8bI5RYzIFDX1jVwyu7DPfdzwONvy4FzAKyZjVNvQQvmkArIyvfjtEltFaQEbG1roCnnyqA0NzRxobe+zWqtHXk4WF88s9OdJEkCso/8eldTS8wJ+QaSKyfXDzHxoFJcOAF4nAAAVs0lEQVQ0FpQWcKKji+0HjoUaR3VdExJcOffMPS4ry4rZvL+V3YePD1Nkri+xzkeSb2YFUa95vau7XN92Hz5BS1unJxKXFHoa2cOem6S6vpHzpoylcHTOGferKo88qFjjgziGKtYSyU2SxkQtj416SNCdQc8fpA+N4pLBnJLRZGdlhNpOcrD1JK/sPkJVWf/jws4qymPa+FHUePVWqGKttL+/5/kOADM7Atwfn5BSS21DC5kZomxiftihONevEZkZlE3ID7VEsmpLE2aRaqv+RLoBF/PC1oO0dSTOg5TpJqZeW/SdcGI9Nq1t2NPMnOLR5I7IDDsU52KyoLSAJ9Y38IllawZ03M3nT+aacyYN+vo19U0U5mVzTozPslSVlbDspZ38bvshrpwX/igW6SjWZLBG0jeIjOZrwJ8Da+MWVQqpbWjh8rlFYYfhXMzeu6iUDQ3N7D58IuZjDh07yXNbmjhv6lhKx579eHJd3cazm5u4qjwy90gs3j6rkJysDKrrGz2RhCTWRPLnwBeBnwbLzwB/E5eIUkjj0TYaj570BxFdUrl0dhFP/fkVAzpm16HjvPMbz/KVX9bxrdvOP+trr991hCPHO2JqH+kxMjuTS2YXUlPfxP3Xn/Wl3SDE2mvrmJnd3TMIopl9wczC7R+YBN4cOt57bLnUNnX8KP73lbN48pUG1uw4dNbnqalvJCOGbr+9VZWVsP3AsdC7LaerWHttrZA0Nmp5nKSn4xdWatgYJJIKTyQuDXyqcjYTC3L50i820n2WDzRW1zdywfRxjBk1YkDH9ZRgvPdWOGLttVUU9NQCwMwO43O296u2oZnphaMoyB3YH4VzyWhUdhZ3X1POa3uaefz3uwd8fGNLGxv2tLxlEqtYTSscxaziPH+eJCSxJpJuSaeGRJE0g9PMA+Le1DMHiXPp4oZFpbxt2lj+4Vf1HG3rGNCxNZsjSWAg7SPRqspKeGnbQU60ezfg4RZrIrkXeF7SjyX9GHiWyLDv7jRa2jrYefC4N7S7tCKJ+69fwIHWk3yneuuAjq2pb2RCQQ7zJ53dM1dVZSW0d3bz0rYDZ3W8O3uxNrb/ClgM1BPpuXUnEHvfwDTk7SMuXZ03dSzvu2AKjzy/gx0xNn53dHXz3OYDVJWVnPUsohfOHMeo7Eyfyz0EsTa2f5zIPCR3Bq8fA38bv7CSX0+PLR8axaWjv353GSMyxQPLN8W0/+93Huboyc6YnmY/nZysTC6dXUR1fSNmXvM+nGKt2voscCGw08yqiEx762n/DGr3NFOSn0Nx/pkHnXMuFZUU5PKZq+awYuN+ntvS/0dFdX0TWRnisjmDe3i3qryY3YdP8HrTGefdc0Ms1kTSZmZtAJJyzKwOKItfWMnPG9pduvvoZTOZNn4US3+xkc6u7jPuW1PfyIUzxpM/yB6OPT2+vHpreMWaSHYHz5H8N7BC0hNAQ/zCSm5tHV1sbWpNiHmvnQtL7ohM7r12PlsaW/nJy2+cdr+GIyeo23f0jJNYxWry2JGUTcj3ya6GWayN7TeZ2REz+1siQ6X8APBh5E+jft9RurrNSyQu7b2rYgKXzSnkGys2c/hYe5/79Dz7cbbdfnurLC9m9Y5DA+5+7M7egOd+NbNnzexJM+v7t8KxIRiC27v+unQnifuuW8DRtg7+6deb+9ynur6RyWNHMqdk9JBcs6qshI4u44WtB4fkfK5/cZ1EXNLVkuolbZV0dx/bp0taKelVSTWSpvTaXiBpj6TvRK2rCc65Pngl3BP2tQ0tFORmMWXc2Y+C6lyqKJuYz+0XT+ffX36D+n1H37LtZGcXL2w9QFV58Vl3++3tgunjyM/J8uFShlHcEomkTCLDzl8DVAC3SarotdvXgGVmdi6wFHiw1/YvE3n4sbfbzWxR8Eq435ZIQ/uYIfvDcC7Z/dUfzWN0ThZLn6p9S9fc1dsPc7y9a8iqtSAyOdcV84qoqW/ybsDDJJ4lkouArWa2LagGewy4odc+FUSeTwGojt4u6QJgApEh65NGZ1c3dXu9x5Zz0cblZfOX75zLC1sPsmLj/lPrq+sbyc7K4JLZhUN6vcqyEva1tFHXqwTk4iOeiWQysCtqeXewLtorwC3B+5uAfEmFkjKArwOfO825Hwmqtb6oBPva/3rTMU52drNgsicS56Ld/vbpzC0ZzQPLN3GyMzIeVnV9I2+fVcio7KGdcLUymODKe28Nj3gmkr4+4HuXM+8ClkhaBywB9gCdwKeB5Wa2iz90u5mdA1wRvD7Q58WlT0paI2lNU9Pw9Smv9YZ25/o0IjOD+66vYOfB4zzywg7eOHicbU3HTn3oD6WSglwWlBZQ48+TDIt4zru+G5gatTyFXs+emFkDcDOApNHALWbWLOkS4ApJnwZGA9mSWoPJtfYExx6V9CiRKrRlvS9uZg8BDwEsXrx42CpKaxtayB2RwayivOG6pHNJ44q5xbxzfgnfXrmFYyc7Aagqj09/maqyEr737Os0H+8Y8PwmbmDiWSJZDcyVNFNSNnAr8GT0DpKKgmosiIwm/DCAmd1uZtPMbAaRUssyM7tbUpakouDYEcB1wIY43sOAbdjTTPnEArIy49ohzrmkde+1FbR3dfPt32xlRuEoZsbpS1dVeTFd3cZzW71UEm9xK5GYWaekO4CngUzgYTOrlbQUWGNmTwKVwIOSDFgFfKaf0+YATwdJJBP4NfD9eN3D15+pZ/2uI/3vGGX9riO874Ip/e/oXJqaWZTHRy+byb+u2nZWk1jFatHUcYwdNYKv/LKOn67uq5a8b1kZ4nPvLveRuwcgnlVbmNlyYHmvdfdFvX8ceLyfc/wQ+GHw/hhwwVDHeTon2rtoDYrfsTpn8hhuPL93nwLnXLQ7rprD7iMnuO2iaf3vfJYyM8QdVXP4n9f2DujvuLahhUkv7+TvbzonbrGlGqVDP+vFixfbmjVrwg7DOZcEPvXjtby6+wgv3H1V2j8LJmmtmS3ubz+vyHfOuShV5cU0NLexeb8PRR8rTyTOORfl1FD0/gxKzDyROOdclAkFuVRMKqC6zhNJrDyROOdcL5VlxazZeZgWH4o+Jp5InHOul6ryErq6jee3HAg7lKTgicQ553o5f+pYCnJ9KPpYeSJxzrlesjIzuHJeMdU+FH1MPJE451wfqspKaDp6ktqGlrBDSXieSJxzrg9LyiKjEnv1Vv88kTjnXB+KRudw3pQxVNf7oI/98UTinHOnUVlWwro3DnP4WHvYoSQ0TyTOOXcaVeUldBus2uKlkjPxROKcc6dx7uQxFOZlU+PVW2fkicQ5504jI0NcOa+YZzc30dXt3YBPxxOJc86dQWVZMYeOtfPq7oFNcpdOPJE459wZXDm3mAzhvbfOwBOJc86dwbi8bM6fNs6fJzkDTyTOOdePqrJiXt3dTNPRk2GHkpA8kTjnXD96Jrtatdmrt/riicQ55/qxoLSAkvwcnzXxNOKaSCRdLale0lZJd/exfbqklZJelVQjaUqv7QWS9kj6TtS6CyS9FpzzW5IUz3twzjlJVJYVs2pzE51d3WGHk3DilkgkZQLfBa4BKoDbJFX02u1rwDIzOxdYCjzYa/uXgWd7rfse8ElgbvC6eohDd865P1BVVkJLWyfrdnk34N7iWSK5CNhqZtvMrB14DLih1z4VwMrgfXX0dkkXABOAZ6LWTQIKzOwli0wSsAy4MX634JxzEZfNLSIrQz6Xex/imUgmA7uilncH66K9AtwSvL8JyJdUKCkD+DrwuT7OubufcwIg6ZOS1kha09TkDWTOucEpyB3BBdPH+fMkfYhnIumr7aL3GAN3AUskrQOWAHuATuDTwHIz29Vr/1jOGVlp9pCZLTazxcXFxQOL3Dnn+lBVXsKmvS3sa24LO5SEEs9EshuYGrU8BWiI3sHMGszsZjM7H7g3WNcMXALcIWkHkXaUD0r6SnDOKWc6p3POxUtV0A3YH058q3gmktXAXEkzJWUDtwJPRu8gqSioxgK4B3gYwMxuN7NpZjaDSKllmZndbWZ7gaOS3h701vog8EQc78E5506ZN2E0pWNyvRtwL3FLJGbWCdwBPA1sAn5mZrWSlkp6b7BbJVAvaTORhvUHYjj1nwH/BmwFXgd+OdSxO+dcXyRRWV7C81sO0N7p3YB7KNL5KbUtXrzY1qxZE3YYzrkUsGLjfj6xbA2PfvxiLp1TFHY4cSVprZkt7m8/f7LdOecG4NLZhWRnZlDjw6Wc4onEOecGIC8ni4tnjffnSaJ4InHOuQGqLCthS2Mruw4dDzuUhOCJxDnnBqiqLPJsmldvRXgicc65AZpZlMe08aOo8eotwBOJc84NmCSqyop54fUDtHV0hR1O6DyROOfcWagsL6Gto5uXtx8KO5TQeSJxzrmzcMmsQnKyMrz3Fp5InHPurOSOyOTS2YU+7haeSJxz7qxVlZew4+Bxth84FnYoofJE4pxzZ6lyXmQ04HSv3vJE4pxzZ2la4ShmF+el/WjAWWEH4JxzyayqrIQfvriD67/9fFyvk5UpPvfuMi6dnXgDRXoicc65Qbjt4mm8ceg4nd3xHUl9Y0MLn/t/r7LyziXkjsiM67UGyhOJc84Nwuzi0Tz0wX5HWh+0F18/wJ98/2W+v2obf/6OuXG/3kB4G4lzziWBS2cXcc3CifxLzesJN2e8JxLnnEsSX3jPfLrM+Oqv6sIO5S08kTjnXJKYOn4Un7xiFv+1bg9rdx4OO5xTPJE451wS+bPK2ZTk57D0F7V0x7mBP1aeSJxzLonk5WRx9zXlvLK7mZ+v2xN2OECcE4mkqyXVS9oq6e4+tk+XtFLSq5JqJE2JWr9W0npJtZI+FXVMTXDO9cGrJJ734JxziebGRZNZNHUs//CrOlpPdoYdTvwSiaRM4LvANUAFcJukil67fQ1YZmbnAkuBB4P1e4FLzWwRcDFwt6TSqONuN7NFwSu9Hyl1zqWdjAxx//UVNB49yb9Ubw07nLiWSC4CtprZNjNrBx4Dbui1TwWwMnhf3bPdzNrN7GSwPifOcTrnXNI5f9o4bn7bZP7tue28cTDcuePj+QE9GdgVtbw7WBftFeCW4P1NQL6kQgBJUyW9Gpzjq2bWEHXcI0G11hclKT7hO+dcYvv81eVkZYoHlm8MNY54JpK+PuB7dzG4C1giaR2wBNgDdAKY2a6gymsO8CFJE4Jjbjezc4ArgtcH+ry49ElJayStaWpqGvzdOOdcgplQkMtnqubwdO1+Xtx6ILQ44plIdgNTo5anANGlCsyswcxuNrPzgXuDdc299wFqiSQNzGxP8PMo8CiRKrQ/YGYPmdliM1tcXFw8NHfknHMJ5mOXz2TKuJEsfWojnV3docQQz0SyGpgraaakbOBW4MnoHSQVSeqJ4R7g4WD9FEkjg/fjgMuAeklZkoqC9SOA64ANcbwH55xLaLkjMvmba+dTt+8o/7F6V/8HxEHcEomZdQJ3AE8Dm4CfmVmtpKWS3hvsVkkkQWwGJgAPBOvnAy9LegV4Fviamb1GpOH96aDtZD2RqrDvx+senHMuGbx7wUQumVXIN56p58jx9mG/vswS48nIeFq8eLGtWbMm7DCccy5uNu1t4dpvPccHL5nB3753wZCcU9JaM+t3aGPvVuuccylg/qQCbrtoGj/+7U627D86rNf2ROKccynir/5oHnnZmSx9aiPDWdvkicQ551JE4egc/uKd83huywF+Uzd8g354InHOuRTygUumM7s4jy8/tZH2zuHpDuyJxDnnUsiIzAzuu34BOw4e54cvbh+Wa3oicc65FLNkXjFXlZfw7ZVbaTp6sv8DBskTiXPOpaC/uXY+F8wYR1tHV9yvlRX3KzjnnBt2s4pH88OP9DmC1JDzEolzzrlB8UTinHNuUDyROOecGxRPJM455wbFE4lzzrlB8UTinHNuUDyROOecGxRPJM455wYlLSa2ktQE7OxntyLgwDCEk4jS+d4hve/f7z19xXL/082suL8TpUUiiYWkNbHMBJaK0vneIb3v3+89Pe8dhvb+vWrLOefcoHgicc45NyieSN70UNgBhCid7x3S+/793tPXkN2/t5E455wbFC+ROOecGxRPJM455wYl7ROJpKsl1UvaKunusOOJN0kPS2qUtCFq3XhJKyRtCX6OCzPGeJE0VVK1pE2SaiV9Nlif8vcvKVfS7yS9Etz7l4L1MyW9HNz7TyVlhx1rPEnKlLRO0lPBclrcv6Qdkl6TtF7SmmDdkP3ep3UikZQJfBe4BqgAbpNUEW5UcfdD4Ope6+4GVprZXGBlsJyKOoE7zWw+8HbgM8H/dzrc/0ngKjM7D1gEXC3p7cBXgW8G934Y+FiIMQ6HzwKbopbT6f6rzGxR1LMjQ/Z7n9aJBLgI2Gpm28ysHXgMuCHkmOLKzFYBh3qtvgH4UfD+R8CNwxrUMDGzvWb2++D9USIfKJNJg/u3iNZgcUTwMuAq4PFgfUreew9JU4BrgX8LlkUa3X8fhuz3Pt0TyWRgV9Ty7mBduplgZnsh8mELlIQcT9xJmgGcD7xMmtx/UK2zHmgEVgCvA0fMrDPYJdV///8J+GugO1guJH3u34BnJK2V9Mlg3ZD93mcNQYDJTH2s8/7QKU7SaOA/gb8ws5bIF9PUZ2ZdwCJJY4H/Aub3tdvwRjU8JF0HNJrZWkmVPav72DUl7x+4zMwaJJUAKyTVDeXJ071EshuYGrU8BWgIKZYw7Zc0CSD42RhyPHEjaQSRJPITM/t5sDpt7h/AzI4ANUTaicZK6vlCmcq//5cB75W0g0gV9lVESihpcf9m1hD8bCTyJeIihvD3Pt0TyWpgbtBzIxu4FXgy5JjC8CTwoeD9h4AnQowlboI68R8Am8zsG1GbUv7+JRUHJREkjQTeSaSNqBp4X7BbSt47gJndY2ZTzGwGkb/z35jZ7aTB/UvKk5Tf8x54F7CBIfy9T/sn2yW9h8g3k0zgYTN7IOSQ4krSfwCVRIaQ3g/cD/w38DNgGvAG8Mdm1rtBPulJuhx4DniNN+vJv0CknSSl71/SuUQaVDOJfIH8mZktlTSLyDf08cA64E/N7GR4kcZfULV1l5ldlw73H9zjfwWLWcCjZvaApEKG6Pc+7ROJc865wUn3qi3nnHOD5InEOefcoHgicc45NyieSJxzzg2KJxLnnHOD4onEuRBImhE9ArNzycwTiXPOuUHxROJcyCTNCubIuDDsWJw7G55InAuRpDIiY399xMxWhx2Pc2cj3Uf/dS5MxUTGN7rFzGrDDsa5s+UlEufC00xkPpzLwg7EucHwEolz4WknMivd05JazezRsANy7mx4InEuRGZ2LJh0aYWkY2aWcsOYu9Tno/8655wbFG8jcc45NyieSJxzzg2KJxLnnHOD4onEOefcoHgicc45NyieSJxzzg2KJxLnnHOD8v8BbrpZ44cNvg0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from week_1_knn import *\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "\n",
    "breast_cancer = load_breast_cancer()\n",
    "\n",
    "X_train = breast_cancer.data[:300, :]\n",
    "y_train = breast_cancer.target[:300, np.newaxis]\n",
    "X_test = breast_cancer.data[300:, :]\n",
    "y_test = breast_cancer.target[300:, np.newaxis]\n",
    "\n",
    "#Plot KNN accuracy against K for K=1-50\n",
    "plot_knn_performance(X_train, y_train, X_test, y_test, maxK=50)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall, the accuracy is very high, so different values for k yield high accuracies. However, arround very high numbers for k (~27) the accuracy starts to drop. \n",
    "In conclusion, for this dataset, the optimum number for k is k=9."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $k$-NN regression\n",
    "\n",
    "Modify the $k$-NN implementation to do regression instead of classification. Compare the performance of the linear regression model and the $k$-NN regression model on the diabetes dataset for different values of $k$.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEWCAYAAABWn/G6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VFX6wPHvm0JCS2iBhBoQBEGKEhAsa/+JFVdRwIqy9rLq2nDtq7urrmV1beyKKCpgN1awIYoKBAydQOiRTiAQICHl/f1xT3SMk2QCmZlk8n6eZ57cuffcc987M5l37r3nniOqijHGGBMqUeEOwBhjTP1iiccYY0xIWeIxxhgTUpZ4jDHGhJQlHmOMMSFliccYY0xIWeIxNUJE2ojIdBHZJSKP10B9o0Tku5qIrT4RkY4iki8i0QGUTRURFZGYCpbfLyKv1XyUkcu99l3CHUdtZ4knAonIahHZJyKtys3PdF80qe55exF5R0S2ikieiCwQkVFuWdmXUn65x/AKNnslsBVIUNW/BG/vfs99QaqInOczL8Z3X2sLF1PXSpaPcmVuKzc/R0SOq6p+VV2rqk1UtaQGwjXV5F77leGOo7azxBO5VgEjy56ISG+gYbkyE4B1QCegJXAJsKlcmWbun6nsMbmC7XUCFut+3JFc0S/uasoFHgzkl/6BEE+w/29ygTtEJCHI26kVAnn/a+gzErT6TPVY4olcE/ASSZlLgVfLlRkAjFfV3aparKo/qeqn1d2QiIx39d/ujopOEpE4EXlKRNa7x1MiEufKH+d+wd8hIhuBlwPYxmMi8p2IJFZQ5DNgH3BRBevHici/RGStiGwSkRdEpKFb1lxEPhKRLSKy3U2391l3mog8LCIzgD1AFxFJFJGXRGSDiPwsIg+VJT0R6Soi37ijyK0iMtnNn+6qnFfF0eMS4Afg5gr2JUpE7hSRFSKyTUTeFJEWbtlvTp+JSGefU6BfiMizfk6fXehel60i8tdyy+JFZLJbf66I9PWJ4xD32uwQkUUicla51+xPPs9/c+rUxXidiCwHlvvZx7L9GC0ia4Gv3PxBIvK92+Y836PAyvZ1P+sbJSIrXX2rRORCN9/v++uzX13ddKKIvOo+V2tE5G5xP1rKXg/3mdzu6j/V3/sdkVTVHhH2AFYDJwFZwCFANL8e2SiQ6sp9AcwARgAdy9WR6srGBLjN8cBDPs8fBH4EWgNJwPfA39yy44Bi4BEgDmjop75RwHd4P47+C0wBGlWw7fuB14CzgJVALBBTbl+fAtKBFkBT4EPgH25ZS+BcoJFb9hbwvk/904C1QC9XbyzwPvAi0Njt4yzgKld+IvBXF3s8cLRPXQp0reR1LNvvfsAOoIWbnwMc56Zvcq9te/f6vQhM9Pe+4SWwfwENgKOBncBr5cr+F+9ouC9QCBzi87oWAcPcPt+KdyQd6x7ZwF2u7hOAXUB3n9fsT+X3q9zr8Ll7P/y9/2Wxvepe44ZAO2AbcJp7bU92z5Oqsa8B1efK7PTZnxSgV3XeX7etD/A+U6nAMmC0z+tRBFyB9/95DbAekHB/f4TkOyrcAdgjCG/qr4nnbuAfwBD3T17+y7g58E9gEVACZAID3LKyf9Qd5R6HVLDN8fw28awATvN5fgqw2k0fh3d0El/JPowCZgKTgXeABpWUvd/nC2am+yf+ZV8BAXYDB/msMxhYVUF9/YDtPs+nAQ/6PG+D9wXd0GfeSOBrN/0qMBZo76fugBKPm34TeMRN+yaeJcCJPuukuC+xGJ/3LQboiJfgG/mUfY3ffxm391k+Cxjh87r+6LMsCtgAHOMeG4Eon+UTgft9XrOqEs8JlbwOZbF18Zl3BzChXLkpeEfbge5roPU1xvu8n0u5xBjI+4uXTAqBnj7LrgKm+bwe2T7LGrl1k0P1PRHOh51qi2wTgAvwPuTlT7OhqttV9U5V7YX3ZZoJvC8i4lOslao283ksCXDbbYE1Ps/XuHlltqhqQRV1dAWGAg+o6r4At3s33q/ReJ95SXj/2HPcKZUdeKfmkgBEpJGIvOhOh+wEpgPN5LfXi9b5THfC+8W/wae+F/GOfABux0t2s9wpqMsDjL28e4FrRCS53PxOwHs+216C98OhTblybYFcVd1TwX6U2egzvQdo4q+8qpbiJcC27rHOzSuzBu8oIlD+YqmsTCfgvLL9dvt+NF7iDXRfA6pPVXcDw4Gr8d7nj0Wkh1svkPe3Fd6RV/n/Ad/X55fX3Sdu39c+YlniiWCqugbv1MhpwLtVlN2Kd5qiLd7pjwO1Hu8fu0xHN++XTQZQxxLgMuBTEekeyEZV9XO8U0DX+szeCuzFO1VSlkATVbXsn/wvQHfgCFVNAP7g5vsmYN941+H9mvVNygkugaOqG1X1ClVti/cr9zmppCVbJfuyFO99u6vconXAqeV+EMSr6s/lym0AWohII595HaoZxi/l3fWJ9njv43qgg/y2oUVHoCyG3XjJvkz55AmBfQbKv+4Tyu13Y1X9J4Hva6D1oapTVPVkvMS2FO+UZKDv71a8o9Dy/wPl36N6yRJP5BuNd0pjd/kFIvKIiBwqXtPjpninqLJVdVsNbHcicLeIJInXrPtevFMf1aKqE/G+eL8QkYMCXO2veL9Ky+ooxfvSeFJEWgOISDsROcUVaYqXmHa4i/T3VRHTBmAq8LiIJIh3sf8gETnW1X2e/No4YTvel11Z8+ZNQHXu83gAL/k285n3AvCwiHRy20sSkaF+4lwDZAD3i0gDERkMnFmNbQP0F5FzxGuscBNewv0R75TmbrwGJbHuovyZwCS3XiZwjjua7Ir3OTxQrwFnisgpIhItIvHiNVRpv5/7WmF94t2XdpaINHb7nI97D6t4fwFQrzn7m3jvU1P3Xt3CfvwPRCJLPBFOVVeoakYFixsB7+Gdy16J9+vsrHJldshv7+O5JcBNP4T3RTAfWADMdfOqTVVfwWus8JUEcF+Oqs7Au1bh6w68I6Ef3em0L/COcsBreNAQ71fqj3in4apyCd6plMV4Xz5v4/0yBq+14EwRycdr0PBnVV3llt0PvOJO7ZwfwL6swjtl2thn9r9dvVNFZJeL+YgKqrgQ73rWNrzXfzLeF2mgPsA75bQduBg4R1WL3KnPs4BT8V6354BL3FEawJN41/E2Aa8Ar1djm36p6jq8U693AVvwjlhu49fvsWrtaxX1ReEdCa/Ha95+LL8eRVf2/vq6AS85r8RrMPIGMK76ex55RLXat10YY+oo1/R3qapWelQXCerTvtY1dsRjTAQTkQHuNGCUiAzB+4X/frjjCob6tK91nd29a0xkS8ZroNASr0XaNar6U3hDCpr6tK91mp1qM8YYE1J2qs0YY0xI2ak2P1q1aqWpqanhDsMYY+qUOXPmbFXVpKrKWeLxIzU1lYyMilogG2OM8UdE1lRdyk61GWOMCTFLPMYYY0LKEo8xxpiQssRjjDEmpCzxGGOMCSlLPMYYY0LKEo8xxpiQssRjjDGGgqISXvhmBXPW5AZ9W3YDqTHG1GMlpco7c3N48vNlbMgr4JrjDqJ/p5oYhLhilniMMaYeUlW+WrqZRz5byrJN+fTt0Iwnzu/H4INaBn3blniMMaaembt2O//8ZCmzVufSuVVjnrvwcE49NBkRCcn2LfEYY0w9sWJLPo99lsVnizbSqkkcD519KMMHdCA2OrSX+y3xGGNMhNu8s4CnvlzO5NnriI+J4paTD2b00Z1pHBeeFGCJxxhjItSugiLGTl/J/75dRXFpKRcP6sT1J3SlVZO4sMZliccYYyJMSanyxsw1PPnFcnJ37+PMvm259f8OplPLxuEODbDEY4wxESVz3Q7ueX8hC37OY3CXltx12iH0bp8Y7rB+wxKPMcZEgB179vHolCwmzlpLUpM4nhl5GGf0SQlZS7XqsMRjjDF1WGmp8vbcHP756VLy9hZx+VGduemkbjSNjw13aBWyxGOMMXXUkg07uef9hWSs2U7/Ts156OxDOSQlIdxhVckSjzHG1DH5hcU8+fkyxn+/msSGsTw6rA/DDm9PVFTtO63mjyUeY4ypI1SVjxds4G8fLWbzrkJGDuzI7ad0p1mjBuEOrVos8RhjTB2wcks+936wiO+yt9KrbQIvXNSfwzo2D3dY+8USjzHG1GLFJaW8OH0l//5iOXGxUTw4tBcXHtGJ6DpyWs0fSzzGGFNLZW/O5y9vzWPeuh2c3juF+87qSeum8eEO64AFtWc4ERkiIlkiki0id/pZHicik93ymSKS6rNsjJufJSKnVFWneB4WkWUiskREbnTzjxORPBHJdI97g7nPxhhzoEpLlZe+W8XpT3/Lmm27eWbkYTx74eERkXQgiEc8IhINPAucDOQAs0UkXVUX+xQbDWxX1a4iMgJ4BBguIj2BEUAvoC3whYgc7NapqM5RQAegh6qWikhrn+18q6pnBGtfjTGmpqzL3cOtb81j5qpcTuzRmn+c2ztiEk6ZYJ5qGwhkq+pKABGZBAwFfBPPUOB+N/028B/xbrMdCkxS1UJglYhku/qopM5rgAtUtRRAVTcHcd+MMaZGqSoTZ63joY8XEy3CY8P6MKx/+1rZ88CBCuaptnbAOp/nOW6e3zKqWgzkAS0rWbeyOg/CO1rKEJFPRaSbT7nBIjLPze/lL1gRudKtm7Fly5bq7KcxxhyQjXkFXPrybO56bwGHdWzGZzf/gfPSOkRk0oHgHvH4e8U0wDIVzfeXKMvqjAMKVDVNRM4BxgHHAHOBTqqaLyKnAe8D3X5XiepYYCxAWlpa+TiNMabGqSrvZ/7MfR8soqhE+ZtrsVZXbgTdX8FMPDl411zKtAfWV1AmR0RigEQgt4p1K5qfA7zjpt8DXgZQ1Z1lhVX1ExF5TkRaqerW/dwvY4w5YFvzC/nrewuYsmgTaZ2a86/z+pLaqnYMWxBswTzVNhvoJiKdRaQBXmOB9HJl0oFL3fQw4CtVVTd/hGv11hnvCGVWFXW+D5zgpo8FlgGISLK7boSIDMTb5201vrfGGBOgzxZu4P+enM7XWVu467QeTL5qcL1JOhDEIx5VLRaR64EpQDQwTlUXiciDQIaqpgMvARNc44FcvESCK/cmXqOBYuA6VS0B8Fen2+Q/gddF5GYgH/iTmz8MuEZEioG9wAiX3IwxJqTyC4t5IH0Rb83JoXe7RB4/vy8Ht2ka7rBCTuw7+PfS0tI0IyMj3GEYYyLI3LXbuXlyJuty93Dd8V258cRuxEYH9VbKkBOROaqaVlU567nAGGOCqLiklOemreDfXy4nOSGeyVcNZkBqi3CHFVaWeIwxJkjW5e7h5smZZKzZztn92vLg2YeSUIsHaAsVSzzGGFPDyppJ3/P+IgT494h+DO1X/jbG+ssSjzHG1KC8vUXc/f5CPpy3ngGpzXni/H50aNEo3GHVKpZ4jDGmhsxcuY1b3pzHxp0F3Pp/B3PNcV3r9PAFwWKJxxhjDlBRSSlPfbGM56atoFOLRrxzzZH069As3GHVWpZ4jDHmAKzYks/NkzOZn5PH8LQO3HtmTxrH2VdrZezVMcaY/bCvuJQXv1nBM19n0zA2mucvPJxTe6eEO6w6wRKPMcZUU8bqXMa8u4Dlm/O9kUHP7EnrhMgaMyeYLPEYY0yA8vYW8ehnS3l95lraNWvIuFFpnNCjTbjDqnMs8RhjTBVUlU8WbOT+DxexLb+Q0Ud35paTD7ZrOfvJXjVjjKnEzzv2cu/7C/ly6WZ6tU1g3KUD6N0+Mdxh1WmWeIwxxo+SUmX896t5fGoWqnD36Ycw6shUYiKsY89wsMRjjDHlLPw5jzHvLmDBz3kc3z2JB4cear0P1CBLPMYY4+zZV8yTny9j3IzVNG/UgP9ccBin907BjSVpaoglHmOMAb5bvpUx781nXe5eRg7syJ1DepDYyHqSDgZLPMaYei1vTxEPfbyYt+bk0KVVYyZfOYgjurQMd1gRzRKPMabe+nTBBu5NX0Tu7n1ce9xB3HhiN+Jjo8MdVsSzxGOMqXc27yzg3g8W8dmijfRqm8DLowZwaDtrIh0qVSYeEUkDjgHaAnuBhcAXqpob5NiMMaZGqSpvZeTw0MeLKSwu5Y4hPbjimM7WRDrEKkw8IjIKuBFYBcwBsoB44GjgDhFZCNyjqmtDEKcxxhyQtdv2MOa9+czI3sYRnVvwz3P70LlV43CHVS9VdsTTGDhKVff6Wygi/YBugCUeY0ytVVKqvDxjFf+amkVsVBQP//FQRg7oSJQN0BY2FR5fquqzFSUdtzxTVb+srHIRGSIiWSKSLSJ3+lkeJyKT3fKZIpLqs2yMm58lIqdUVad4HhaRZSKyRERu9Jn/tCs/X0QOryxmY0zkWLpxJ+c8N4OHPl7C0V1bMfWWP3DhEZ0s6YRZtRoXiMhcVQ3oi1tEooFngZOBHGC2iKSr6mKfYqOB7araVURGAI8Aw0WkJzAC6IV3bekLETnYrVNRnaOADkAPVS0Vkdau/Kl4R2bdgCOA591fY0yE2l1YzNNfLuel71aR2DCWZ0Yexhl97EbQ2qK6rdqq864NBLJVdSWAiEwChgK+iWcocL+bfhv4j3ifjKHAJFUtBFaJSLarj0rqvAa4QFVLAVR1s882XlVVBX4UkWYikqKqG6qxL8aYOkBV+XThRv720WI25BUwPK0Dd5zagxaNG4Q7NOOjuonn42qUbQes83mew++PNH4po6rFIpIHtHTzfyy3bjs3XVGdB+EdLf0R2ALcqKrLK4ijHWCJx5gIsnJLPvelL+Lb5VvpmZLAfy44nP6dmoc7LONHZa3axB0l/EJV766qjO9iP/PKl62oTEXz/V2TKqszDihQ1TQROQcYh9cMPJA4EJErgSsBOnbs6GcVY0xttHdfCc9+nc3Y6SuJi4ni/jN7ctGgTtZEuhar7IjnaxF5B/jAt8m0iDTAa1J9KfA1ML6C9XPwrrmUaQ+sr6BMjojEAIlAbhXrVjQ/B3jHTb8HvFyNOFDVscBYgLS0tIqSqTGmFvl88SYe+HAROdv38sfD2jHmtB60bmpDUNd2lf0kGAKUABNFZL2ILBaRVcByYCTwpKqOr2T92UA3EensktUIIL1cmXS8BAYwDPjKHUGlAyNcq7fOeA0DZlVR5/vACW76WGCZzzYuca3bBgF5dn3HmLptXe4eRo+fzRWvZtCoQTSTrhzEk8P7WdKpIyo84lHVAuA54DkRiQVaAXtVdUcgFbtrNtcDU4BoYJyqLhKRB4EMVU0HXgImuMYDuXiJBFfuTbxGA8XAdapaAuCvTrfJfwKvi8jNQD7wJzf/E+A0IBvYA1wWSPzGmNqnoKiEsdNX8uzX2URHCXed1oPLjupMrJ1Wq1Ok4ks09VdaWppmZGSEOwxjjI9vlm3hvg8WsnrbHk7vk8Ldpx9CSmLDcIdlfIjIHFVNq6qcdRJqjKnV1uXu4W8fLWbq4k10btWYCaMHcky3pHCHZQ6AJR5jTK1UUFTC89NW8MI3K4iOEu4Y0oPLj04lLsaGLajrKk08rveBKap6UojiMcbUc6rK1MWb+NtHi8nZvpcz+7blrtN62Gm1CFJp4lHVEhHZIyKJqpoXqqCMMfXTii35PPDhYqYv20L3Nk2ZeMUgBh9ko4FGmkBOtRUAC0Tkc2B32UxVvTFoURlj6pX8wmKe+Wo5475bRXxsNPed2ZOL7SbQiBVI4vmY6nWVY4wxAVFV0uet5++fLGHTzkLO69+eO07tQasmceEOzQRRlYlHVV9xN2uW9Q6dpapFwQ3LGBPplmzYyX3pi5i1Kpfe7RJ5/qL+HN7R+larDwIZ+vo44BVgNV6/Zx1E5FJVnR7c0IwxkWh3YTFPfL6Ml2d4Qxb845zenJ/WgWgbI6feCORU2+PA/6lqFoAbF2ci0D+YgRljIs/nizdx3wcLWZ9XwAVHdOT2U7rTrJENWVDfBJJ4YsuSDoCqLnNd6BhjTEA25hVwX/pCpizaRPc2TXnHhiyo1wJJPBki8hIwwT2/EJgTvJCMMZGipFSZ8MNq/jV1GUUlpdw+pDtXHNPF+lar5wJJPNcA1wE34l3jmY7XeagpZ9mmXXw0bz2jj+5CYiM7KDT126L1edz17gLm5eRxTLdWPHx2bzq2bBTusEwtEEjPBS+p6kXAE6EJqe5atXU3T3+Vzf/1SiaxUWK4wzEmLHYXFvPUF8sYN2M1zRvF8vTIwzizTwreqPbGBNZzQZKINFDVfaEKqq5KSfTGAtmQV8Ch7SzxmPrnq6WbuOf9Rfy8Yy8jB3bgziGH2NG/+Z1ATrWtBmaISDq/7bnAjoDKSU7wEs/GvL1hjsSY0Nq0s4AHPlzEJws20q11E966ejADUluEOyxTSwWSeNa7RxTQNLjh1G0tm8QREyVsyCsIdyjGhERpqfL6rLU8+ulSCktKue0Ur/FAgxhrPGAqFsg1niaqeluI4qnToqOENgnxbLTEY+qBZZt2MebdBcxZs52jurbk4bN7k9qqcbjDMnVAINd4Dg9VMJEgOTHejnhMRCsoKuG5r7N5/psVNImL4fHz+nLO4e2s8YAJWCCn2jLd9Z23+O01nneDFlUdlpwYz5L1O8MdhjFB8ePKbdz13gJWbtnNOYe146+nH0JL69DTVFMgiacFsA04wWeeApZ4/EhJiOerJZtRVfsFaCJG3p4i/vHpEibNXkeHFg159fKB/OFgG37a7J9Aeqe+LBSBRIrkxHj2FpWwc2+xNSM1dZ6q8tH8DTzw4WK279nHVcd24aYTD6ZhAxt+2uy/CpueiMibPtOPlFs2NZhB1WXJZffy7LQm1aZuy9m+h9GvZHDDxJ9o2yye9OuPYsyph1jSMQessiOebj7TJwN3+Dy3Y+wK+N5E2iM5IczRGFN9JaXK+O9X8/hUr2/ge87oyagjU23YAlNjKmtsr/u57BciMkREskQkW0Tu9LM8TkQmu+UzRSTVZ9kYNz9LRE6pqk4RGS8iq0Qk0z36ufnHiUiez/x7A4l9fyUnNgSwJtWmTtq8s4DzX/yBv320mCM6t2DqzX9g9NGdLemYGlXZEU8jETkMLzk1dNPiHg2rqtjdA/Qs3tFSDjBbRNJVdbFPsdHAdlXtKiIjgEeA4SLSExgB9ALaAl+4cYCoos7bVPVtP+F8q6pnVBVzTWjdNA4RrEm1qXMy1+3gqgkZ7NxbzFPD+zG0X1trIGOCorLEs4FfOwbdyG87Cd0YQN0DgWxVXQkgIpOAoYBv4hkK3O+m3wb+I94nfSgwSVULgVUiku3qI4A6wyo2OoqkJnFsssRj6pB35uQw5r0FtG4ax7vXHskhKXaa2ARPhYlHVY8/wLrbAet8nucAR1RURlWLRSQPaOnm/1hu3XZuurI6H3an0r4E7nSJC2CwiMzD6/rnVlVdtN97FYCUxHg27LTEY2q/4pJS/v7JUsbNWMXgLi159sLDadHYRgQ1wRXMDpX8HaOXvzZUUZnqzgcYA/QABuDde1TWGGIu0ElV+wLPAO/7DVbkShHJEJGMLVu2+CsSsOTEeOso1NR623fv49KXZzFuxipGHZnKq6MHWtIxIRHMxJMDdPB53h7viMNvGRGJARKB3ErWrbBOVd2gnkLgZdypOVXdqar5bvoTIFZEWpUPVlXHqmqaqqYlJR1Yo73kBOs2x9RuWRt3MfTZGcxetZ1Hh/Xh/rN62aigJmSC+UmbDXQTkc4i0gCvsUB6uTLpwKVuehjwlaqqmz/CtXrrjNe0e1ZldYpIivsrwNnAQvc82c1DRAbi7fO2IO0z4LVs21VQTH5hcTA3Y8x++WzhBv743AwKikqYdNUgzk/rUPVKxtSgCq/xVNU5qKrOrWJ5sYhcD0wBooFxqrpIRB4EMlQ1HXgJmOAaD+TiJRJcuTfxGg0UA9epaomL63d1uk2+LiJJeKfjMoGr3fxhwDUiUgzsBUa45BY0ZffybMwroGvrJsHclDEBKy1VnvpyOU9/uZx+HZrx4sX9aePGkDImlKSi72AR+dpNxgNpwDy8L/U+wExVPTokEYZBWlqaZmRk7Pf6P67cxoixP/La6CM4utvvzuoZE3L5hcXcPDmTzxdvYlj/9jx09qHEx1oPBKZmicgcVU2rqlyVrdpck+UrVXWBe34ocGtNBRqJfjnisZZtphZYvXU3V7yawcqtu7nvTK8XArs/x4RTIL1T9yhLOgCqurCsVwDjXxsbAtvUEp8v3sRf3swkKkp49fKBHNXVjsBN+AWSeJaIyP+A1/CaLl8ELAlqVHVcfGw0LRo3sJZtJmzy9hTxwEeLeHfuz/RMSeCFi/rTsWWjcIdlDBBY4rkMuAb4s3s+HXg+aBFFCBsC24TLtKzN3PHOfLbm7+PGE7py/QndaBBjTaVN7RHIeDwFIvIC8ImqZoUgpoiQYkNgmxDbVVDEwx97g7V1a92E/16SRp/2zcIdljG/U2XiEZGzgMeABkBnd33nQVU9K9jB1WXJifFkrtsR7jBMPTEjeyu3vz2fDXl7ufrYg7jppG7Was3UWoGcarsPrxeAaQCqmuk7fIHxLyUhntzd+ygoKrEvABM0uwuL+eenS5nw4xq6tGrMW1cfSf9OzcMdljGVCiTxFKtqnjW/rJ6ykUg37yy0i7omKGatyuXWt+axbvseLj+qM7ed0t1GBzV1QiCJZ6GIXABEi0g34Ebg++CGVfeluAHhNuTttcRjalRBUQmPTcli3IxVdGjeiElXDOKILi3DHZYxAQsk8dwA/BUoBN7A667moWAGFQmS7SZSEwRz127n1jfnsXLrbi4e1Ik7T+1B47hA/o2NqT0q/cS6UUQfUNXb8JKPCVBZ4rGWbaYmqCovz1jNQx8vJiWxoXXHZOq0ShOPqpaISP9QBRNJmsTF0DQuxu7lMQdsX3Ep936wkEmz13FyzzY8fn5fEuJjwx2WMfstkGP0n0QkHXgL2F02U1XfDVpUESI5MZ4N1m2OOQDb8gu55rW5zFqdy3XHH8RfTu5OVJQ19DF1WyCJpwXe+DUn+MxTwBJPFZIT49m4s7Dqgsb4sXTjTkaPz2BrfiH/HtGPof3aVb2SMXVAID0XXBaKQCJRSmI8yzYd2DDapn6aumgjN0/OpHFcDG9eNZi+HawHAhM5Aum5IB4YDfQbDQD8AAAgAElEQVTCG5sHAFW9PIhxRYTkxIZs3lVIUUmpDStsAqKqPDdtBf+amkWfdomMvSTNBmszESeQb8MJQDJwCvAN0B7YFcygIkVyQjyqsGWXnW4zVSsoKuGmyZk8NiWLM/u0ZfJVgy3pmIgUSOLpqqr3ALtV9RXgdKB3cMOKDCnWpNoEaNPOAoa/+AMfZK7ntlO68+8R/ayrJROxAmlcUOT+7nCjj24EUoMWUQT55SZSSzymEvPW7eDKCRnsKijmxYv7c0qv5HCHZExQBZJ4xopIc+AeIB1oAtwb1KgihA2BbaqSPm89t701j1ZN4njnmiM5JCUh3CEZE3SBtGr7n5v8BugS3HAiS2LDWOJjo2wIbPM7xSWlPPH5Mp6btoIBqc154aL+tGwSF+6wjAmJQFq1+T26UdUHaz6cyCIipCQ2tGs85jc27Szgxok/MXNVLiMGdODBoYfaCKGmXgnkVNtun+l44AxgSXDCiTzJNgS28fHNsi3cMjmTPftKePy8vpzbv324QzIm5Kr8maWqj/s8HgaOAwK6hVpEhohIlohki8idfpbHichkt3ym7wBzIjLGzc8SkVOqqlNExovIKhHJdI9+br6IyNOu/HwROTyQ2GtKsg2BbfBOrT02ZSmXjptFyyYN+PCGoyzpmHprf/pTb0QA13pcz9bPAicDOcBsEUlX1cU+xUYD21W1q4iMAB4BhotIT2AE3k2rbYEvRORgt05ldd6mqm+XC+VUoJt7HAE87/6GRHJiPJt2FlBaqtbHVj21Mc87tTZrdS7D0zpw/1m9bMA2U68Fco1nAV7fbADRQBIQyPWdgUC2qq509UwChgK+iWcocL+bfhv4j3hDnQ4FJqlqIbBKRLJdfQRQZ3lDgVdVVYEfRaSZiKSo6oYA9uGApSTGU1yqbN1dSOumdjNgfTMtazO3vDmPgqISnhzelz8eZkc5xgRyxHOGz3QxsElViwNYrx2wzud5Dr8/0viljKoWi0ge0NLN/7HcumWn9yqr82HXGOJL4E6XuPzF0Q74TeIRkSuBKwE6duwYwO4FJtndeb4pzxJPfVJcUsrjny/j+Wkr6JHclP9ccDhdWzcJd1jG1AqBNKXZ5fPYCySISIuyRyXr+TuvpAGWqe58gDFAD2AAXo/ad1QjDlR1rKqmqWpaUlKSn1X2j+8Q2KZ+2JC3l5H//ZHnp61g5MAOvH/dUZZ0jPERyBHPXKADsB3vS7wZsNYtUyq+3pPj1ivTHlhfQZkcEYkBEoHcKtb1O9/n1FmhiLwM3FqNOILGhsCuX75euplb3sxkX3GpDWVgTAUCOeL5DDhTVVupaku8U2/vqmpnVa2skcFsoJuIdBaRBniNBdLLlUkHLnXTw4Cv3LWYdGCEa/XWGa9hwKzK6hSRFPdXgLOBhT7buMS1bhsE5IXq+g5Ay8YNiI0Wa9kW4YpKSvnHp0u4bPxskhMb8uENR1vSMaYCgRzxDFDVq8ueqOqnIvK3qlZy12yuB6bgNUoYp6qLRORBIENV04GXgAmu8UAuXiLBlXsTr9FAMXCdqpYA+KvTbfJ1EUnCOyrLBMpi/gQ4DcgG9gAhHV8oKkpo3dTu5YlkOwuKuOa1OczI3sYFR3Tk3jN6WgefxlRCvAOMSgqITAG+BV7DO7V2EfAHVT2l0hXrsLS0NM3IyKix+oY9/z0x0cKkKwfXWJ2mdli/Yy+XvTybFVvy+cc5vTkvrUPVKxkToURkjqqmVVUukFNtI/GaUL8HvA+0dvNMgLx7eWxMnkizeP1O/vjcDNbv2Mv4ywZa0jEmQIF0EpoL/BnA9VK9Q6s6TDK/kZIYzxdLNqGqeJegTF03fdkWrn19Lk3jY3jrmsH0SLZepY0JVIVHPCJyr4j0cNNxIvIV3nWSTSJyUqgCjATJiQ0pKColb29R1YVNrffm7HVcNn427Zs35L1rj7KkY0w1VXaqbTiQ5aYvdWVbA8cCfw9yXBHFRiKNDKrKE1OzuP2d+Rx5UEveunrwL83ljTGBqyzx7PM5pXYKMFFVS1R1CfvXx1u9ZSOR1n37ikv5y1vzePqrbM5Pa8+4UQNoGh8b7rCMqZMqSyCFbqjrTcDx/HpDJngdhZoAlXWbY0c8ddPOgiKunjCH71ds45aTD+aGE7ratTpjDkBliefPeB13JgFPquoqABE5DfgpBLFFjKSmcUQJNhJpHeTbXNrGzzGmZlSYeFR1Jl7fZ+Xnf4J3U6YJUGx0FElN46zbnDpm0fo8Lh8/mz2FJbxy+UCO6toq3CEZExHsWk2IJNsQ2HXKN8u2cO1rc0hoGGvNpY2pYZZ4QiQlIZ4VW/LDHYapgqry0ner+MenSzm4TVNeHjXAWq4ZU8Ms8YRIcmI8M7K3hjsMU4mdBUXc9tY8pizaxCm92vCv8/payzVjgiCgxCMiRwKpvuVV9dUgxRSRUhLj2VVYzK6CIvsyq4UWr9/Jta/PYd32vdx9+iGMPrqztVwzJkgCGfp6AnAQXo/PJW62ApZ4qqHsdM2mnQWWeGqZtzLWcff7C0lsGMukKwcxILWy8Q2NMQcqkCOeNKCn9c92YHzv5enaummYozEABUUl3PfBIiZnrGNwl5Y8PfIwkprGhTssYyJeIIlnIZAMhGzwtEhUNgS29V5QO6zZtptrXpvL4g07ue74g7jl5O5ER9mpNWNCIZDE0wpYLCKzgF/69lfVs4IWVQRqneD9krbEE35TF23kL2/NI0qEcaPSOKFHm3CHZEy9EkjiuT/YQdQH8bHRtGzcgA12E2nYFJeU8tiULF6cvpLe7RJ57sLD6dDCen8yJtQCGY/nm1AEUh8kJ9oQ2OGyeWcB10/8iVmrcrnwiI7cY8NTGxM2gbRqGwQ8AxwCNACigd2qardyV1NKYjw/77DEE2o/rtzG9W/8xO7CYp44vy/nHG79rRkTToEMff0fvKGulwMNgT+5eaaa2iTEW0ehIfbx/A1c/NJMEhrG8P51R1nSMaYWCOgGUlXNFpFoVS0BXhaR74McV0RKSYxn+54iCopK7DRPCEyatZa73lvA4R2b89KoASQ2tPunjKkNAkk8e0SkAZApIo/iNatuHNywIlOya1K9aWcBnVraSxhMY6ev4O+fLOXYg5N44aL+NGxgid6Y2iKQU20Xu3LXA7uBDsC5wQwqUtkQ2MGnqjz62VL+/slSTu+Twn8vSbOkY0wtU2XiUdU1gAApqvqAqt6iqtmBVC4iQ0QkS0SyReROP8vjRGSyWz5TRFJ9lo1x87NE5JRq1PmMiOT7PB8lIltEJNM9/hRI7MFgQ2AHV2mpcs8HC3lu2gpGDuzA0yMOo0FMIL+tjDGhVOV/pYiciddP22fueT8RSQ9gvWjgWeBUoCcwUkR6lis2Gtiuql2BJ4FH3Lo9gRFAL2AI8JyIRFdVp4ikAc38hDNZVfu5x/+qij1YbAjs4CkqKeXmNzN57ce1XHVsF/7+x97WE4ExtVQgPwfvBwYCOwBUNROvp+qqDASyVXWlqu4DJgFDy5UZCrzipt8GThSvS+ChwCRVLXRDbme7+iqs0yWlx4DbA4gtLBrHxZAQH2Mt22pYQVEJV02YwweZ67l9SHfGnHqI9SxtTC0WSOIpVtW8/ai7HbDO53mOm+e3jKoWA3lAy0rWrazO64F0VfXXp9y5IjJfRN4WkQ7+ghWRK0UkQ0QytmzZEsj+7ZfkxHg74qlBuwqKuGTcLL7O2sxDZx/Ktcd1DXdIxpgqBJJ4ForIBUC0iHQTkWeAQJpT+/vJWb6H64rKVGu+iLQFzsO70bW8D4FUVe0DfMGvR1i/rUR1rKqmqWpaUlKSvyI1IjmxIZus25wasS2/kJH//ZG5a7bz1PB+XDSoU7hDMsYEIJDEcwPetZZCYCKwE7gpgPVy8FrAlWkPrK+ojIjEAIlAbiXrVjT/MKArkC0iq4FGIpINoKrbVLWsc9P/Av0DiD1oUhLsiKcmrN+xl/Nf/IHlm/IZe0l/hvYrfzBtjKmtAumrbQ/wV/eojtlANxHpDPyM11jggnJl0oFLgR+AYcBXqqqu8cIbIvIE0BboBszCO+L5XZ2qughv6AYARCTfNVhARFJ8Tr+dBSyp5n7UqOTEeLbkF1JUUkpstLW42h+rtu7mov/NZOfeIiaMPoKBnW3gNmPqkgoTT1Ut16oaFkFVi0XkemAKXv9u41R1kYg8CGSoajrwEjDBHZ3k4iUSXLk3gcVAMXCd6zUBf3VWsY83ishZrp5cYFQV5YMqJTEeVdi8q5B2zRqGM5Q6aeHPeYx6eRaqMPHKQRzaLjHcIRljqkkqGlhURLbgXcifCMyk3PWVSO61Oi0tTTMyMoJS97SszYx6eTbvXDOY/p3sl3p1TJ69lns+WESrxg2Y8KcjOCipSbhDMsb4EJE5qppWVbnKTrUlAyfjdRB6AfAxMDGAIwxTibKRSO06T+AKikq494OFvJmRw9FdW/HvEf1o2cSGqDamrqow8bhTW58Bn4lIHF4CmiYiD6qqv9ZjJgBlN5Fa7wWB8R2i+sYTuvLnkw62G0ONqeMqbVzgEs7peEknFXgaeDf4YUWuhIYxNIyNtsQTgCmLNnLrW/OIjhJevmwAx3dvHe6QjDE1oLLGBa8AhwKfAg+o6sKQRRXBRISUxHgbArsSxSWlPDY1ixe/WUmf9ok8e4ENUW1MJKnsiOdivN6oD8ZrGVY2XwC1EUj3nw2BXbHNuwq44Y2fmOmGqL73zJ7ExVjv0sZEksqu8dhNJkGSnBjPzJW54Q6j1pm5chvXT/yJXQVFPDm8L388zEYLNSYSBTQCqalZKYnxbNpZQEmp2oVyvDF0xk5fyaNTsujUohGvjT6C7slNwx2WMSZILPGEQXJiQ4pLlW35hbR2rdzqq7y9Rdz21jymLt7Eab2TeeTcPjSNtyGqjYlklnjC4Jcm1TsL6nXiWbZpF1e8msHP2/dyzxk9ufyoVBvOwJh6wBJPGPgOgd2nnl7G+HrpZm6Y+BMNG0Qz6cpBpKVaLw7G1BeWeMKgPg+Braq89N0q/v7JEg5JSeB/l6b90puDMaZ+sMQTBi0aNaBBdFS96zZnX3Ep936wkEmz1zGkVzJPDO9Lowb2ETSmvrH/+jCIihLaJMbVqyGwc3fv45rX5jBzVS7XH9+VW04+mChr0WdMvWSJJ0xSEhrWmyOe5Zt2MfqVDDbuLODfI/rZoG3G1HOWeMIkOTGeeTk7wh1G0E3L2swNb/xEXGw0k68cxGEdm4c7JGNMmFnvBGFS1m1OReMh1XWqyrjvVnH5+Nm0b9GID64/ypKOMQawI56wSU6Ip7C4lB17imjeuEG4w6lRRSWl3PvBIibOWsv/9WzDk8P70TjOPmrGGI99G4SJ7708kZR4tu/ex7Wvz+WHldu47viD+MvJ3a0RgTHmNyzxhMkv9/Ls3EvPtpHR0Xf25nxGvzKbDTsKeOL8vpxzeD29O9YYUylLPGESaUNgZ2/exbAXfiAmSph45SD6d7LrOcYY/yzxhElS0ziioyQiei9Yv2MvF780i9joKN6+ejCdWjYOd0jGmFrMWrWFSXSU0LppXJ1PPDv27OPScbPILyjmlcsGWtIxxlTJjnjCqE1CPBvr8BDYe/eVcPn42azJ3cOrlw+MmGtVxpjgCuoRj4gMEZEsEckWkTv9LI8Tkclu+UwRSfVZNsbNzxKRU6pR5zMikh/INsItJTG+zl7jKSop5drX55C5bgdPj+jHoC4twx2SMaaOCFriEZFo4FngVKAnMFJEepYrNhrYrqpdgSeBR9y6PYERQC9gCPCciERXVaeIpAHNAtlGbVB2E2ldo6rc8c58vs7awkNn92bIoSnhDskYU4cE84hnIJCtqitVdR8wCRharsxQ4BU3/TZwongjgQ0FJqlqoaquArJdfRXW6ZLSY8DtAW4j7FIS48kvLGZXQVG4Q6mWf366lHfn/swtJx/MBUd0DHc4xpg6JpiJpx2wzud5jpvnt4yqFgN5QMtK1q2szuuBdFXdEOA2fkNErhSRDBHJ2LJlS4C7eGCSXZPqunTU89/pK3lx+kouGdyJG07oGu5wjDF1UDATj7+jivIdk1VUplrzRaQtcB7wzH7GgaqOVdU0VU1LSkrys0rN8+29oC54d24OD3+yhNN7p3Dfmb1smGpjzH4JZuLJATr4PG8PrK+ojIjEAIlAbiXrVjT/MKArkC0iq4FGIpJdxTbCLjmhrPeC2p94vs7azO1vz+eori15Ynhfoq0bHGPMfgpm4pkNdBORziLSAK+xQHq5MunApW56GPCVet01pwMjXIu0zkA3YFZFdarqx6qarKqpqpoK7HGNCSrbRti1SagbQ2DPXbuda1+bS4+UprxwUX/iYqLDHZIxpg4L2n08qlosItcDU4BoYJyqLhKRB4EMVU0HXgImuKOTXLxEgiv3JrAYKAauU9USAH91VhGK323UBg1iomjVpEGtPtWWvXkXl4+fTeuEOF4eNZCm8bHhDskYU8dJLfnxX6ukpaVpRkZGSLZ1xjPfktQkjpcvGxiS7VXH+h17Gfb89xSVKu9cfSQdWzYKd0jGmFpMROaoalpV5azLnDBLTmhI1sZd5BcWhzuU3yjrCmdXQTHjLxtgSccYU2Ms8YTZhYM6smlXIZe/PJs9+2pH8lm+aRfnPv89a3L3MPaSNHq1TQx3SMaYCGKJJ8yO796ap4b3I2NNLle8mkFBUUlY4/kg82eGPjuDvL1FvHLZQAYfZF3hGGNqliWeWuDMvm15bFhfvl+xjWtem0NhceiTT2FxCfd+sJA/T8qkV9sEPr7xGEs6xpigsN6pa4lz+7dnX0kpY95dwA1v/MSzFx5ObHRofhf8vGMv174+l3nrdnDFMZ25fUiPkG3bGFP/WOKpRUYO7Mi+4lLuS1/ETZMz+ffwfsQEOQF8s2wLN036iaIS5YWLDrcOP40xQWeJp5a59MhU9hWX8vAnS4iLjuJf5/UlKgi9BJSUKk9/uZynv1pO9zZNef6i/nRuZYO4GWOCzxJPLXTFH7pQWFzCv6Yuo0FMFH//Y+8aTT65u/fx50k/8e3yrZxzeDsePrs3DRtYbwTGmNCwxFNLXX9CNwqLS3nmq2waxETxwFk10ynnT2u3c93rc9m6ex//OKc3IwZ0sM4+jTEhZYmnFrvl5IMpLC5l7PSVNIiO4q+nH7LfSUJVmfDjGv720WLaJMTzztVH0ru93Z9jjAk9Szy1mIgw5tQeFBaV8L/vVhEfG82tp3SvVh35hcUs+jmP12au5cN56zmhR2ueOL8vzRo1CFLUxhhTOUs8tZyIcN+ZvdhXUsp/vs4mLiaKG07s5rfs3n0lLN6wkwU5O5j/cx7zc/JYsSUfVYgSuO2U7lxz7EFBaaxgjDGBssRTB0RFCQ+f3ZvColIe/9xrcDDqqFSyNu5ifk4e83N2MD8nj+Wb8ykp9Tp9TWoaR9/2iZzZpy192ifSp30iLZvEhXlPjDHGEk+dERUlPDqsD/tKSvnHp0v519Qsikq8JNO8USy92zfj5J5t6N0ukT7tm9EmIc4aDRhjaiVLPHVITHQUTw7vR9fWTdhbVEKfds3o0z6R9s0bWpIxxtQZlnjqmNjoKG466eBwh2GMMfvNOuQyxhgTUpZ4jDHGhJQlHmOMMSFliccYY0xIWeIxxhgTUpZ4jDHGhJQlHmOMMSFliccYY0xIiaqGO4ZaR0S2AGuqKNYK2BqCcGqj+rzvUL/3vz7vO9Tv/Q9k3zupalJVFVni2U8ikqGqaeGOIxzq875D/d7/+rzvUL/3vyb33U61GWOMCSlLPMYYY0LKEs/+GxvuAMKoPu871O/9r8/7DvV7/2ts3+0ajzHGmJCyIx5jjDEhZYnHGGNMSFni2Q8iMkREskQkW0TuDHc8wSQi40Rks4gs9JnXQkQ+F5Hl7m/zcMYYLCLSQUS+FpElIrJIRP7s5teX/Y8XkVkiMs/t/wNufmcRmen2f7KINAh3rMEiItEi8pOIfOSe16d9Xy0iC0QkU0Qy3Lwa+exb4qkmEYkGngVOBXoCI0WkZ3ijCqrxwJBy8+4EvlTVbsCX7nkkKgb+oqqHAIOA69x7XV/2vxA4QVX7Av2AISIyCHgEeNLt/3ZgdBhjDLY/A0t8ntenfQc4XlX7+dy/UyOffUs81TcQyFbVlaq6D5gEDA1zTEGjqtOB3HKzhwKvuOlXgLNDGlSIqOoGVZ3rpnfhfQG1o/7sv6pqvnsa6x4KnAC87eZH7P6LSHvgdOB/7rlQT/a9EjXy2bfEU33tgHU+z3PcvPqkjapuAO/LGWgd5niCTkRSgcOAmdSj/XenmjKBzcDnwApgh6oWuyKR/Pl/CrgdKHXPW1J/9h28HxlTRWSOiFzp5tXIZz+mhgKsT8TPPGuTHsFEpAnwDnCTqu70fvjWD6paAvQTkWbAe8Ah/oqFNqrgE5EzgM2qOkdEjiub7adoxO27j6NUdb2ItAY+F5GlNVWxHfFUXw7Qwed5e2B9mGIJl00ikgLg/m4OczxBIyKxeEnndVV9182uN/tfRlV3ANPwrnU1E5GyH62R+vk/CjhLRFbjnU4/Ae8IqD7sOwCqut793Yz3o2MgNfTZt8RTfbOBbq51SwNgBJAe5phCLR241E1fCnwQxliCxp3TfwlYoqpP+CyqL/uf5I50EJGGwEl417m+Boa5YhG5/6o6RlXbq2oq3v/4V6p6IfVg3wFEpLGINC2bBv4PWEgNffat54L9ICKn4f36iQbGqerDYQ4paERkInAcXpfom4D7gPeBN4GOwFrgPFUt3wChzhORo4FvgQX8ep7/LrzrPPVh//vgXUCOxvuR+qaqPigiXfCOAloAPwEXqWph+CINLneq7VZVPaO+7Lvbz/fc0xjgDVV9WERaUgOffUs8xhhjQspOtRljjAkpSzzGGGNCyhKPMcaYkLLEY4wxJqQs8RhjjAkpSzxmv4lIqm+v1ZWUS/Hp3fc4EVEROdNn+Uc+d4cfaEyrRaRVTdRVxXYecz02PxbsbflsM01Enq5gWVD2O9D3OMC6HhSRk6ooc7+I3BrMOAIhIg1EZLrPzaKmBtmLakLhFuC/Ps9zgL8CH4YnHP9EJManH66qXAUkhfIeDlXNADJqoi4RiXbd4YSMqt4byu0FoqL3XFX3iciXwHDg9dBHFtnsiMfUCBHp4sYtGeBn8bnAZz7P5wF5InKyn3p++eXufuFPc9P3i8grIjLVlTlHRB5144V85rq2KXObG0dmloh0desnicg7IjLbPY7yqXesiEwFXi0Xi7gjm4VuO8Pd/HSgMTCzbJ7POk1E5GVXfr6InOvmj3TzForIIz7l80XkYfHGvPlRRNq4+ee5svNEZLqbd5zPkWNL91r8JCIv4tOPmIhc5PY9U0ReFG8oj7JtPSgiM4HBItJfRL4RrxPIKfJrVyj93XZ/AK7z836WxTJNRN4WkaUi8rqI14ldJfWOF5Fhbvo0t953IvJ02X45PV3dK0XkRp/5Me4zMN9tt5Gr60T3OiwQb/youAA+S7+85yLSy+f1mi8i3dz23gcu9Lf/5gCpqj3ssV8PIBWvG43ueHdx9/NTpjMwx+f5ccBHwDHAN27eR8Bxbno10MpNpwHT3PT9wHd4XfP3BfYAp7pl7wFn+6z/Vzd9CfCRm34DONpNd8TrBqes3jlAQz+xn4vXI3M00AbvTu0Utyy/gtfkEeApn+fNgbZu3SS8swxf+cSrwJlu+lHgbje9AGjnppv5vnZu+mngXjd9uqunFV4nnh8CsW7Zc8AlPts6303HAt/jHbWB98t+nJueDxzrph8DFvrZz+OAPLz+yqKAH4Cjq6h3PF53M/F4Pbx3dvMn+uzX/W79OLc/21ydqS7+o1y5ccCtPnUd7Oa/iteZa9lnoaLP0i/vOfAMcKGbbuAzPxrYEu7/s0h82BGPOVBJeP01XaSqmX6WpwBbys9U1W8BROSYamzrU1UtwvtSjubXo6gFeF9MZSb6/B3spk8C/iNeF//pQIK4vqiAdFXd62d7RwMTVbVEVTcB3wD+juh8nYQ3UCAAqrrdrTNNVbeod1rndeAPrsg+vMQL3pdh2X7MAMaLyBVuX8v7A/Ca28bHeIOSAZwI9Admu309EejilpXgdXgK3o+FQ/F6Hc4E7gbai0giXqL7xpWbUMm+zlLVHFUtBTJd7H7rLbdeD2Clqq5yzyeWW/6xqhaq6la8TijbuPnrVHWGm34N7/3pDqxS1WVu/iv8+tpWxvc9/wG4S0TuADqVzVfvVOQ+n8+JqSF2jcccqDy8X5xHAYv8LN+L96vUn4fxrvX4nmMv5tdTwOXXKwRQ1VIRKVL3sxSvHzXfz7L6mY4CBpdPMO7s0O4K4tuf8Q+E33eVX1k9vvtRgtsPVb1aRI7AO5rJFJF+ftb119+VAK+o6hg/ywr01+s6AixS1cG+BcTrFDTQfrR8r2+Vxe63Xj8xVrde/MSlVdRV2Wfpl/dcVd9wpx9PB6aIyJ9U9Su3OA4oqCJeU012xGMO1D68UQgvEZEL/Cxfxm+PRn6hqlPxTkX19Zm9Gu8XO3inuvbHcJ+/P7jpqcD1ZQUq+CIvbzowXLzB0JLwfknPqmKd8ttpjtep6LEi0spdbxmJd/RUIRE5SFVnqndBfiu/HYqjLLYLXdlT8V5H8IYjHibeGCqISAsR6eRnE1lAkogMduViRaSXesMf5InXQSpU/xqH33rLlVkKdBFvcD349f2qSseyevFew+9cXaniruUBF/Pra7uaAD5L4nWIuVJVn8Y7Gu7j5rfEO9VWFGB8JkCWeMwBU9XdwBnAzSIy1M+yFT5fDOU9zG9PxTwA/FtEvsX7tbs/4twv2D8DN7t5NwJp7uLxYuDqAOp5D+96xzy86zK3q+rGKtZ5CGhe1jAAb8z6DcAYvC715wFzVbWq7uQfcxfLF+IlmXnllj8A/EFE5uJ1Wb8WQFUX453emioi8/GuUaWUr+TyLswAAADcSURBVFy9YduHAY+4ODOBI93iy4BnXeMCf6cgK1RFvWVl9gLXAp+JyHd4vZ7nBVD9EuBSt18tgOdVtcDF+5aIlPUi/oIrH+hnaTiw0J0a7MGvjUyOBz4JIC5TTdY7tQk6Efkj0F9V7w53LKZ2EJEmqprvWsI9CyxX1SfDHZcvEXkXGKOqWeGOJdLYEY8JOlV9D++0hzFlrnBHGIuARP6/vTsmAgCEgSCY+FeACySmoEIAT7OrIN3NV6lan++59HnyuEXnDYsHgCiLB4Ao4QEgSngAiBIeAKKEB4CoAeKdSqabSsC3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from week_1_knn import *\n",
    "\n",
    "regression_plot() #activate the function to apply kNN regression and plot resulting MSE's"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the selected values for k, which are all odd numbers from 1 until 50, show a significantly lower Mean Squared Error (MSE) compared to the linear regression model used above. The minimum MSE is achieved by considering 5 neighbours. The entire graph indicates that the kNN regression model is significantly better at performing the task at hand compared to the linear regression model. This can be explained, since the considered neighbours are all selected by taking the correct feature values into account. The nearest neighbours are hence close to the original test point in the feature space, while the points further away are not considered. Taking the mean of these neighbours, which are known to be close to the true test point, does hence provide a similar result. This algorithm is therefore biased, as the original feature values have a significant importance and influence in producing the final prediction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class-conditional probability\n",
    "\n",
    "Compute and visualize the class-conditional probability (conditional probability where the class label is the conditional variable, i.e. $P(X = x \\mid Y = y)$ for all features in the breast cancer dataset. Assume a Gaussian distribution.\n",
    "\n",
    "<p><font color='#770a0a'>Based on visual analysis of the plots, which individual feature can best discriminate between the two classes? Motivate your answer.</font></p>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from scipy.stats import norm\n",
    "%matplotlib qt\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "\n",
    "breast_cancer = load_breast_cancer()\n",
    "\n",
    "#Extract data\n",
    "X = breast_cancer.data\n",
    "Y = breast_cancer.target[:, np.newaxis]\n",
    "\n",
    "#Normalize data to [0,1] range\n",
    "X_norm = (X-X.min(axis = 0))/(X.max(axis = 0)-X.min(axis = 0))\n",
    "\n",
    "#Set boundaries for the gaussian distributions and create subplots\n",
    "x = np.linspace(norm.ppf(0.05), norm.ppf(0.95), 100)\n",
    "fig, axs = plt.subplots(5, 6, constrained_layout=True)\n",
    "\n",
    "def gaussian(x, mu, sig):\n",
    "    return (1/(sig*(2*math.pi)**.5))*np.exp((-(x - mu)**2.)/(2*sig**2.))\n",
    "\n",
    "#Draw subplots for each feature\n",
    "for feat in range(X.shape[1]):\n",
    "    feat_data = X_norm[:,np.newaxis,feat]\n",
    "    #Extract data with corresponding label 0\n",
    "    axs.flat[feat].plot(x, gaussian(x,feat_data[Y==0].mean(),feat_data[Y==0].std()),'r-', lw=2, alpha=0.6, label='norm pdf')\n",
    "    #Extract data with corresponding label 1\n",
    "    axs.flat[feat].plot(x, gaussian(x,feat_data[Y==1].mean(),feat_data[Y==1].std()),'b-', lw=2, alpha=0.6, label='norm pdf')\n",
    "    axs.flat[feat].set_title(str(feat+1))\n",
    "    axs.flat[feat].set_xlabel(\"x\")\n",
    "    axs.flat[feat].set_ylabel(\"P\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How well a feature can seperate the classes can be derived from the degree in which the gaussian distributions shown in the figure are seperated. It appears that the 28th feature shows the least overlap between the distributions (relative to their sizes) and can therefore best discriminate between the two classes. "
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}