3. model_selection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "=================================================================================================================\n",
    "# Lecture Notes: Model Selection\n",
    "\n",
    "\n",
    "##### D.Vidotto, Data Mining: JBI030 2019/2020\n",
    "\n",
    "\n",
    "=================================================================================================================\n",
    "\n",
    "As you may have noticed, in the past lecture we have trained *k-nearest-neighbors* and *logistic regression* models on a subset of the dataset called *training data*, and evaluated it with several performance measures on the *test data*. In this lecture, we are going to explore why this step is necessary. In particular, we are going to introduce a very important concept in Data Mining and statistical learning in general: the *Bias-Variance* tradeoff. Understanding this concept is crucial to understand how to perform model selection, and how to tune what are called the *hyperparameters* of the models we will explore in the rest of the course. We will see two strategies that will help us to perform the tuning step: the *Holdout method* and *Cross-Validation*.  \n",
    "<br>\n",
    "This lecture is structured as follows:\n",
    "\n",
    " 1. some important definitions \n",
    "   * signal/noise \n",
    "   * overfitting/underfitting\n",
    "   * bias-variance tradeoff\n",
    "   * model complexity\n",
    "   * model hyperparameters\n",
    "   * training/testing errors\n",
    " 1. model selection and validation strategies \n",
    "   * Holdout method \n",
    "     * standard split\n",
    "     * stratified split\n",
    "   * Cross-Validation\n",
    "     * standard cross-validation\n",
    "     * stratified cross-validation\n",
    "     * leave-one-out cross-validation \n",
    "   * grid search and random search  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Some Important Definitions \n",
    "\n",
    "**Signal and Noise**. \n",
    "\n",
    "<br>\n",
    "\"<i>The signal is the truth. The noise is what distracts us from the truth.\"</i><sup>1</sup><br>\n",
    "<sup>1</sup>N.Silver (2012). <i>The Signal and the Noise.</i> Penguin Group.\n",
    "\n",
    "<br>\n",
    "\n",
    "When we observe the raw data, we must be aware that not all the information contained in such data is relevant. Indeed, we can decompose the data in this way: \n",
    "\n",
    "$$Data = Relevant + Irrelevant Information$$\n",
    "\n",
    "which can also be rewritten as \n",
    "\n",
    "$$Data = Signal + Noise$$\n",
    "\n",
    "Therefore, we can define the **Signal** as that part of relevant information that is present in the data (the one that is of interest for our business goal/research question), and that we want to uncover. Conversely, the **Noise** is the (unavoidable) component of irrelevant information (for our goal at hand) that is present in the data. In order to discover the signal, we implement different Data Mining models on the dataset, hoping that one of them will result useful for our purposes.\n",
    "***\n",
    "**Overfitting and Underfitting**. However, when we train an algorithm on a dataset, it is not certain at all whether what we are fitting is the right model; and the same can be said even when we try multiple algorithms. In particular, there are two types of mistakes that an algorithm can do when training a model: \n",
    " \n",
    " 1. it can ignore important aspect of the data, that are relevant for our business goal; \n",
    " 1. it can learn too closely from the given data points, so that it doesn't know how to generalize to new ones. \n",
    " \n",
    "In general, in case (1) the noise is disregarded by the model, but the model is so \"concerned\" to ignore it that even (part of) the signal is discarded. This is also known as **model bias**, and in general it happens when we are trying to train a too simple model for the data and research question at hand. On the other hand, in case (2) the model is too concerned to learn so many aspects of the data, that it learns the data points 'by heart'; thus, when presented -in the model implementation stage- to new data, it fails to correctly recognize how to behave with them, which increases the likelihood to perform wrong predictions. When the model is discarding relevant part of the signal in the data is said be **underfitting**; when, on the other hand, the model is learning the training examples too closely is said to be **overfitting**. Of course, there exist models in the middle between under- and over-fitting ones; those models should be the ones achieving best performance when performing predictions on a new set of data, and are the ones that we are seeking. \n",
    "***\n",
    "**Bias-Variance tradeoff**. In case (1) above, the model will produce in general wrong (but rather stable across different samples) predictions (hence it is a *low variance, high bias* case). On the other hand, models in case (2) will fit very closely the training data (hence low bias), but the predictions become highly unstable (high variance) if trained on a different sample, and won't generalize well on new data. We are going to see an example of biased and unstable algorithms shortly.\n",
    "\n",
    "Ideally, we would like to minimize both bias and variance simultaneously. However, this is not possible, as models with decreasing[increasing] bias have increasing[decreasing] variance. This can be explained mathematically. The *expected MSE* (mean squared error) of a new point $(x^*,y^*)$, predicted by an already trained model to be $\\hat{f}(x^*)$ when in reality it should be $f(x)$, is (the solution requires a couple of simple algebraic steps): \n",
    "\n",
    "<br>\n",
    "\n",
    "$$ \\mathbb{E}[(y^* - \\hat{f}(x^*))^2] = (Bias[\\hat{f}(x^*)])^2 + Var[\\hat{f}(x^*)] + Var(e^*)$$\n",
    "\n",
    "\n",
    "where:\n",
    "\n",
    " * $Bias[\\hat{f}(x^*)] = \\mathbb{E}[\\hat{f}(x^*)]-\\mathbb{E}[f(x^*)]$ is the model **bias** introduced above; it captures that part of error given by the simplifying assumptions of the model w.r.t. the true function (for instance, when modelling non-linear data (true function:$f(x^*)$) with linear models (estimated function:$\\hat{f}(x^*)$) (Bias is part of the *reducible error*; it can be lowered by using a more complex model or perhaps by increasing the number of features)\n",
    " \n",
    " \n",
    " \n",
    " * $Var[\\hat{f}(x^*)] = \\mathbb{E}\\{(\\hat{f}(x^*)-\\mathbb{E}[\\hat{f}(x^*)])^2\\}$ represents the **variance** of the estimator itself; that is, how much would the prediction vary across different training samples? (For example, linear models are in general rather stable, while more complex models have higher variance) (Variance is part of the *reducible error*; it can be lowered by adding more training points or by reducing model complexity)\n",
    " \n",
    " \n",
    " \n",
    " * $Var(e^*) = \\mathbb{E}\\{(y-\\mathbb{E}[f(x^*)])^2\\}$ is the *irreducible error*, which represents the **noise** (=the difference between the true function $f(x^*)$ and the actual observation); this is that part of the error that cannot be captured by the model, under which the error test cannot lie, as it is intrinsic in the data  (If this error is large, it might be lowered only with data cleaning, adding new features to the dataset, and so on).\n",
    " \n",
    " \n",
    "This formula is also known as *Bias-Variance decomposition* of the expected error. It shows us how the expected test error is decomposed, and it represents a *tradeoff* between model accuracy (bias) and model stability (variance). This tradeoff is known as **Bias-Variance Tradeoff**. Models with high variance (and low bias) are said to be **overfitting**, while models with high bias (low variance) are usually **underfitting**. The goal in Data Mining is to find a model that is in between, that is, it is not too biased, but at the same time not too unstable (as mentioned above, unfortunately it is not possible to obtain both simultaneously). Consider the following examples in which we try to predict a continuous output $y$ with a numeric feature $x$ (thus we are in a *regression setting*): \n",
    "\n",
    "\n",
    "\n",
    "<img src=\"./img/model_selection/bias_variance_1.png\" width=\"800\" />\n",
    "\n",
    "\n",
    "There are three models trained (fitted) to the same dataset: a linear regression model, a polynomial regression model with degree 3 (a *polynomial regression* model simply exponentiate the feature $x$ to a given degree, and then fit a linear model on such new features); and a polynomial regression model with degree 10. Could you recognize which model is underfitting, which one is overfitting, and which one represents the best fit? <br>\n",
    "\n",
    "\n",
    "The \"true\" signal in the preceding example was a polynomial of degree 3 (since data were artificially simulated, the true generating function was known). It is clear, then, that a linear regression appears to make too strict assumptions on the dataset, while a polynomial of degree 10 is left too \"free\" to vary. \n",
    "\n",
    "\n",
    "We are now going to see how *biased models* (left plot above) remain stable across different samples, while *unstable models* (right plot above) can vary significantly their predictions across different samples. The following plots show two different datasets (samples) generated from the same signal as before, on which the same underfitting and overfitting models used above are trained:   \n",
    "<img src=\"./img/model_selection/bias_variance_2.png\" width=\"700\" height=\"50\"/>\n",
    "\n",
    "\n",
    "As you can see, the simple model is consistenlty biased (it doesn't pick up all relevant aspects in the signal), but it always fits a straight line (although with different slopes and intercept); the model on the right, instead, is more wiggly, and it tries to learn every time  different aspects of the dataset, resulting in dramatically different predictions. What we are observing here in a regression setting also happens in classification: the simpler the classifier, the \"smoother\" the decision boundary will be (and the higher the probability of underfitting); the more complex the classifier, the more wiggly and close to the data points the decision boundary. \n",
    "\n",
    "Whether a model is under/over-fitting is also a function of the sample size: the larger the training sample size, the more complex the model is \"allowed\" to be, and the true function will be recovered correctly (however, the higher the complexity of the algorithm we train, the larger the sample size required to retrieve the true $f(x)$). In the following plots, you can see the same data generated in the previous example, where a polynomial regression of degree 10 is fitted with different sample sizes: n=35, n=100 and n=150.  \n",
    "<img src=\"./img/model_selection/bias_variance_3.png\" width=\"700\" height=\"1000\"/>\n",
    "\n",
    "Therefore, the larger the sample size, the more an overfitting model can \"adjust\" towards the true function. In reality, however, we don't know the true generating function, nor we know what is the correct sample size for the model we are fitting; therefore, we need to operate in a different manner. The way to go is by \"manually\" regulating the *model complexity*. \n",
    "\n",
    "***\n",
    "\n",
    "**Model Complexity**. So far we have encountered the term *complexity* of a model, without really defining it. From the plots seen above, it should be quite intuitive now. Model complexity refers to the *degree of flexibility* of a model, that is, how much the model is able to learn. Simplest models try to memorize very basic aspects of the data (such as a linear function, or a linear decision boundary). More complex models, instead, are capable of memorizing more complex relationships. As we have seen above, however, too complex models run into the risk of learning too closely from the training data points (without generalizing well to new ones). Therefore, increasing model complexity corresponds to increasing its variance (and reducing bias), while decreasing model complexity corresponds to increasing its bias (and reducing variance). \n",
    "\n",
    "In general, increasing the complexity of the model also means decreasing its interpretability, and vice versa. *Interpretability* of a model refers to the human ability to interpret the results of a Machine Learning algorithm; for example, interpreting the regression coefficients of a linear regression model. For this reason, highly complex models (which could be most suitable for some datasets) are also known to be *black-box* models, as they are of difficult interpretation; for such algorithms, we will mostly be interested in their black-box predictions.\n",
    "\n",
    "But how do we determine the complexity of a model? It turns out that each model can be controlled by the experimenter by regulating its *hyperparameters*. \n",
    "\n",
    "*** \n",
    "\n",
    "**Model Hyperparameters**. The hyperparameters of a model are parameters that control model complexity. They are called like this because they operate *on top* of the normal model parameters; that is, the specific values of the hyperparameters of a model will affect the values of the trained parameters, and consequently the model predictions (for *non-parametric* models, hyperparameters directly influence the predictions). In these notebooks and during the lectures, sometimes we may call the hyperparameters with the name of *parameters*, as a short-hand; we will generally refer to hyperparameter, though. The distintiction will be clear from the context.\n",
    "\n",
    "The polynomial degree chosen for the input features above is an example of hyperparameter; as you could have noticed, modifying such hyperparameter drastically changed the predictions of the model. Another example of model hyperparameter is the number of neighbors $k$ in k-Nearest-Neighbors (a model introduced during the last lecture, and that we will encounter again later on). A further example of hyperparameter is the number of input features used to train a regression model. \n",
    "\n",
    "Let's see now an example of how the hyperparameters affect the predictions in a classification example. Once again, we will use synthetic data to show this principle. The example below show the different *decision regions* and *decision boundaries* estimated by k-Nearest-Neighbors with varying number of neighbors: 1, 15, and 50. \n",
    "<img src=\"./img/model_selection/knn_dec_bound.png\" width=\"1500\" height=\"100\"/>\n",
    "\n",
    "Lower values of $k$ lead to a more \"wiggly\" decision boundary, while larger values lead to a smoother one. What do you think that this means? Does increasing $k$ lead to a simpler or more complex models? \n",
    "\n",
    "As we will see, each algorithm has its own set of hyperparameters. Hyperparameters of a model is not something that is learned by the algorithm in the *training stage* (see next subsection), but it is something that we, as users, must set beforehand. Our role as data scientists is to try out different values of the hyperparameters, and find out which one leads to the best performance: such optimal value, of course, can vary from dataset to dataset. To be able to tune the hyperparameters we need to be able to understand the difference between *training* and *testing* stage. \n",
    "\n",
    "***\n",
    "**Training and Test Errors**. In the previous lecture, we have used a *training dataset* to train the algorithms (training means learning the parameters of the algorithm, given a value of the hyperparameters), and a separate *test set* to evaluate how it generalizes (for example, with the RMSE or the $f_1$ score). This is a crucial point in Data Mining and Machine Learning: if we evaluated the algorithm with the same set of data points used for training (called *training set*), we would be able to detect when the algorithm is underfitting, but not when the algorithm is underfitting. This is because, as we increase model complexity, the error on the training set (which, by the way, is known as **training error**) becomes smaller and smaller, approaching $0$ (or, to be more precise, the irreducible error). This makes sense, intuitively: an overfitting model learns so closely the training data points that their error keeps decreasing. On the other hand, an underfitting model never fits well the data points, leading to a large error also when evaluated on the training data . \n",
    "\n",
    "In contrast, assessing the performance of the algorithm with a new set of data (called *test set*) can help us to identify when the algorithm is generalizing well. This step is called *testing stage* or *inference*. In particular, given a specific value of the hyperparameters, we can understand whether we are in an under/over-fitting regime by comparing the performance on the test set (which provides us with an estimate of the **test error**) and training set. \n",
    "\n",
    "The following graph represents the relationship between model complexity, training error, test error, bias and variance: \n",
    "\n",
    "<img src=\"./img/model_selection/test_train_errors.png\" width=\"700\" height=\"100\"/>\n",
    "\n",
    "From the plot we can infer the following points: \n",
    "\n",
    "* when the algorithm is underfitting, the training error is close to the test error, and they are both large; in this case, model bias is large and model variance is low. The model doesn't learn enough from the training dataset, and therefore cannot generalize well either\n",
    "\n",
    "* when the algorithm is overfitting, the training error is small and the test error is large again; here, variance is small and bias is large. The model learns too closely the training data, and cannot generalize well to the test points \n",
    "\n",
    "* the test error has an approximate U-shape: it is large at the extremes, where either bias or variance are large (i.e., when the model does not learn enough or when it learns too much from the training data), but it reaches a minimum (a \"sweet spot\") somewhere in the middle, when the levels of bias and variance both reach a reasonable compromise; ideally, this is the model we should be able to detect, and that can be achieved with a proper tuning of the hyperparameters\n",
    "\n",
    "* the errors leading to this shape of the curve are measures those that we seek to minimize: miss-classifications, RMSE, MAE, and so on. For measures that we want to maximize, such as accuracy, $R^2$, $f_1$, etc. the plot will be symmetrical ($\\cap$-shape rather than $U$-shape of the metric curve on the test set), but the principles explained here still hold\n",
    "\n",
    "* the training error should (at least in theory) always be larger than the test error; this makes sense, as the model model learns better the data points used in the training stage than a set of points that has never seen before \n",
    "\n",
    "In the next section we are going to see how to perform model tuning (a.k.a. model selection) and validation. Let's conclude here by observing how training and test error behave on a simulated dataset. In particular, we are going to implement it on a two-classes example, with 1000 examples (80% of which will be retained in the training set, and the reamining ones will go to the test set) and 30 features (only 3 of such features are actually significant predictors for the classification task). We will try values for $k$ ranging between 1 and 100. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dd3xUVfbAv4dQggQBAUUpC4LrbhBEYFXAgtixK1iAn4pi7wVFRVHUXQXXsooFsSEKoujKKnbFgiJFI1IFFAVFmnRDSXJ+f9z38maSmckkmckkmfP9fO5n3rvv3vvOK3PPu+eWI6qKYRiGkb7USLUAhmEYRmoxRWAYhpHmmCIwDMNIc0wRGIZhpDmmCAzDMNIcUwSGYRhpjimCJCIiy0TkqFTLUV5EZJGIHFoJ5BgnInemWo6qioj8XUTWxjieJSIqIs0qUq54EJHPROSMJJR7oogsTHS5VQ1TBNUMEZknIlu8kC8i20L2by1Lmaq6r6p+nmhZE4WI3B5yjdu86/b3vytHuSVWEiLymohsDznfFhH5sqznTCaqukBVm/j7IjJLRM5OpUzxoqqHqeqkijyn9/zfqshzpgpTBNUMVW2vqlmqmgV8Dlzp76vqP4umF5GaFS9lYlHVu0Ou+Urg85Br3r8CRBgWcr4sVe0eKVGke13a+18dnldRquM1VTVMEVQQXrP8JxE5x9tfJiI3isgcEdkoIq+ISKZ3rKeIrBCRG0RktYisFJGBCZJjkNfM/o+I/AEMFZF9ROQTEflDRNaKyIsi0iAkzwoR6elt3yMi4z0zzWYRmSsinWOc7zEv/yYRmSki3UOOxSxLRLqISI53bDxQpxzX3VFEporIehGZLyInhRw73TN/bRaR5SJyhYjsDrwK/DXkS79B9DNEPOd+XgvlUhFZAUyOFOelPVNEFnjyfSAibUPKWSsi14vIfOCPCOf5t4j8y9vOEpGdIjLM22/knW8X/9xe/MPAAcDz3rXdH1LkiSLyo/c+PBDl2tp5+eqFxB3qPesaIpLtvWd/eO/wcyKSFe2aROQuEXmhyDmeFZF7ve3C1ouIXCki74vIKO+/s0REjgjJt6+IfOU9zykiMkZEnozzmd3ivXN7FImvKSJPiMga75w5ItIunjKrBKpqIUkBWAYcBXQGfgFOLHJsBrAXsBuwALjUO9YTyAOGA7WA3sCfQCPveD9gThznnwoMKhI3yCv7MiADqAv8FTgSqA3sDkwDHgjJswLo6W3fA+QCx3r5RwJfxJDh/7zrqwncDPwK1CmpLFylvwK42rsHZwM7gTtLuOZBwNQicQ2B370yMoCDcRVqG0CADUAXL20ToJO3fSKwsITzvQYMiXJsP6AAeMq7z3WjxHUCNgOHec/gLuB7IMMrZy0wHdgTqBvhPCcDX3nbxwBLgY9Cjk0LkWdbSL5ZwNkh+1mA4hRgfaAtsAk4JMr1zQDOCdl/wn9vgGzce1zLk3sGcE9I2rBrCjlXXe94XWAjsG9RWXGtvp24/0EGcCOwxDsmwBzvHtbGvdd/Ak9GuYbCZwzcD3yF9z8rku4MXAt7V9wHdAegaarrmESFlAtQnQOusr+LkIq0yLEBIfsj/JfV+wPlAjVDjq8GDi7l+acSWRH8WEK+PsDMkP2iiuDdkGMdgS1xyiO4Cq99SWUBvYDlgIQcn0HZFMGFwDtF4l4CbvBkWgucC2QVSROvIsjFKRM/PO4d2w9Xse4ekj5S3P3AsyH7NXGKqqu3vxY4M4YMjYBtQD3vnl6HU3y1gX8D/ww5dzyKoFNI3BSceTHSea8FJofIvAboHCXtAJzJjmjXhKtoz/G2zwZmRJIVpwhyQo7t7smdhVNAW4DaIcf/S2xF8APwJPAhUC9KupNxyvkfoe9kdQlmGko+lwJfqurUCMd+D9n+E/ci+6xT1bwYx8vD8tAdEWkmIhNF5FcR2QQ8j/syjkZRuetFSygiN4nIQhHZCKz30oaWHa2svYAV6v0LPX6OIVMs/gL0EpENfgBOAfb0yj8FOAtYLiIfxTJ1ReEuVW0YEi4PObZDVVcXSV80bi9Crs177r8BzUPShD2zUFR1PbAQ6IFrVXwMfIOrtA4HPi3l9cR6L0N5BThaRBrhWnVrVfUbABFpIa4j/TfvnXqS4u9U0Wt6Aacw8H5fLIWMeHLuBaxW1R0xzlOUZsB5wN2qujVKmv958j0NrPJMnruUUG6VwRRB8rkUaCUiD6VakBCKLjl7P7Ad6KCquwLn476Uy4Vnt70e16xuiPty3RJn2SuBFkXiWpVRlOXAlCKVdZaq3gigqtNU9QRgD1wl+pKXLxFL80Yqo2jcbzhlBRR2nu6FM6PFKieUT3GVcTvcl+unuK/Y/XCmvnhlixtVXQl8iXu+/YCXQw7/G2fayfbeqUsp/tyLnn8icKiItMe1CCeUQayVwO4iUiskrmUJeX4DTgcmikiXSAnU8YCqdgL2B7oCV5VBvkqJKYLksxk4DjhMRO5LtTBRqA9sBTaKSEuczTVR5ebhzAC1gDuJ0XoowhdADa9jsKaInInraykLrwMHikgfr6zaItLN6/CsLyJniUh9nN15C86GD7AK2CO0QzRJTAD6iEgPrwK7DVehfVuKMj4FLsKZ9ApwZsHLgO9UdUuUPKuAvcssteNlnOntZMIVQX3cu79JRFrjzFUxUdVNwNvAeOBDVV1TBnnmAz8Bt4pILe9j5Og4zv0O7v69LSLFRpqJSHdxgxdq4t6RHQTvSZXHFEEFoKobcC/j8SJyd3nLE5H+IjKv/JIVMgw4EPcFNxlI1HjtKTi762Jcn8gmXAVXIqq6HTgN9+dc723/tyxCqOo63NfyRTiTwm8EHfHg+hWW4+z7Z+NaRACzgXeBXzyTUrRRQ3dJ+DyCZaWU71vgEuAZnJ39MOBUVc0vRTGf4Uwjn4XILiH7kfg3cKF3bf8qjcwhTMKNPpqvqktD4ocCR+Ce+Wu4Duh4eAHXERvLLBQVz9R3JnA87r0Z7J1/exx5J+P6H971WiWhNALG4t6RH3Ed8o+VRcbKiISbYA3DMFKHiGTjzE3NVHVbgsp8GzeAYGQiyquOWIvAMIxKgYhk4ExIY8ujBDyzXysRyRCRU3FDuCcnSs7qiM3oMwwj5YibwPcTzox4bDmLa4kzBzXEzd85T1UXlbPMao2ZhgzDMNIcMw0ZhmGkOVXONNSkSRNt3bp1qsUwDMOoUsyePXutqjaNdKzKKYLWrVsza9asVIthGIZRpRCRqDPzzTRkGIaR5pgiMAzDSHNMERiGYaQ5Va6PwDCMimXnzp2sWLGCbdsSMtHXSDKZmZm0aNGCWrVqlZzYwxSBYRgxWbFiBfXr16d169aIlHtRWiOJqCrr1q1jxYoVtGnTJu58ZhoyDCMm27Zto3HjxqYEqgAiQuPGjUvdejNFYBhGiZgSqDqU5VmlpSKYNQv69YM1ZVnt3DAMo5qRlopg8GAYPx6efz7VkhiGURLr1q2jU6dOdOrUiWbNmtG8efPC/R07dpRcADBw4EAWLYq97tyoUaN46aWXYqaJl0MOOYR99923UM6zzjorIeUmi7TrLP79d/jU8+A6L5GuXQzDSAqNGzcmJycHgDvvvJOsrCxuvDHciV6hE/Yakb9tn3vuuRLPc8UVV5Rf2BBeeeUVOnXqFPV4Xl4eNWvWjLofb75EkHYtgtdeA3/B1blzUyuLYRhlZ8mSJWRnZ9O/f3/at2/PypUrufjii+natSvt27dn+PDhhWkPOeQQcnJyyMvLo2HDhgwZMoT999+fbt26sXr1agCGDh3Kww8/XJh+yJAhHHjggey77758+eWXAGzdupUzzjiD7Oxs+vTpQ9euXQuVVDwMGDCAyy67jAMPPJBbb72VoUOHcu6559KjRw/OP/98cnNzOe+88+jQoQOdO3fms8+cg7kxY8Zw6qmncsQRR3DsseVdpbs4adcimDgx2F6wAAoKIMpHhGEYRUlWp3EZl8NfuHAhY8eOpWvXrgDcd9997LbbbuTl5XHEEUfQp08fsrOzw/Js3LiRww8/nPvuu4/rr7+eZ599liFDhkQQSZkxYwaTJ09m+PDhvPvuuzz66KM0a9aMSZMm8d1339G5c3Q32meddRZ169YF4LjjjuO++5zL8pUrVzJ9+nRq1KjB0KFDWbhwIZ999hmZmZncf//91KlTh++//5558+bRu3dvFi9eDMC3335LTk4OjRo1KtO9ikVaKYJff4UvvoDMTKhXD9atg2XLYO/yuu82DCMltG3btlAJAIwfP55nnnmGvLw8fvvtN+bPn19MEdStW5fjjz8egC5duvD5559HLPv0008vTLNs2TIAvvjiC26++WYA9t9/f9q3L+raOCCaaahv375hJqxTTjmFzMzMwvIHDx4MQPv27dlrr71YsmQJAMccc0xSlACkmWnINwv17g0HHODirJ/AMEqBanJCGalXr17h9uLFi3nkkUf4+OOPmTNnDscdd1zE8fS1a9cu3M7IyCAvLy9i2XXq1CkxTXlljrQfb75EklaKwDcLnXkm+Irc+gkMo3qwadMm6tevz6677srKlSt57733En6OHj16MNGrSL7//nvmz5+f0PIPPfTQwpFLCxYsYOXKlbRr1y6h54hE2piGli+HL7+EunXhhBNg0yYXby0Cw6gedO7cmezsbP72t7/xl7/8hR49eiT8HFdddRXnnnsu2dnZhaFBgwYR04b2Eeyxxx5xKaarrrqKSy65hA4dOlCrVi3Gjh0b1oJJFlXOZ3HXrl21LI5pHnwQbrgB+vZ1LYOvvoLu3aFTJ/j22yQIahjVhAULFvD3v/891WJUCvLy8sjLyyMzM5PFixdzzDHHsHjx4oQP5ywvkZ6ZiMxW1a6R0lcu6ZPIK6+43zPPdL9+/9HChZCfDxkZqZHLMIyqw5YtWzjyyCPJy8tDVXnqqacqnRIoC1X/CuLkoYdcZ3Hv3m6/QQNo0QJWrIAff4R99kmtfIZhVH4aNmzI7NmzUy1GwklqZ7GIHCcii0RkiYgUH6gbpDtDRFREIjZbEkH37s48tMsuQZx1GBuGYSRREYhIBjAKOB7IBs4RkewI6eoD1wBfJ0uWaOy3n/u1DmPDMNKZZLYIDgSWqOqPqroDmACcEiHd3cD9QIW7P/JbBKYIDMNIZ5KpCJoDy0P2V3hxhYhIZ6Clqr4dqyARuVhEZonIrDUJXDvaFIFhGEYKO4tFpAbwIHB+SWlVdTQwGtzw0UTJEDpyaOdOKIWLT8MwKoh169Zx5JFHAvD777+TkZFB06ZNAZgxY0bc4+yfffZZevfuTbNmzYodGzBgANOmTSucE1C/fv2oS09UR5KpCH4FWobst/DifOoD+wFTPY86zYDJInKyqpZ+okAZyMqC1q3dekNLloANlTaMykc8y1DHw7PPPkvnzp0jKgKAhx56iFNPPTVq/sq0bHSiSaZ0M4F9RKQNTgGcDfTzD6rqRqCJvy8iU4EbK0oJ+LRv7xTB6NFu/aFatdzM4113rUgpDMMoCy+88AKjRo1ix44ddO/enccee4yCggIGDhxITk4OqsrFF1/MHnvsQU5OTuFs33hbEkOHDuWXX35h6dKltGnThp49e/LWW2+xceNGatSowQcffMCNN97I+++/j4gwbNgw+vTpw4cffsg999xDVlYWS5cuZcGCBRVwN8pO0hSBquaJyJXAe0AG8KyqzhOR4cAsVZ2crHOXho4d4e23wVuGHID+/WHcuNTJZBiVlcq0CvXcuXN54403+PLLL6lZsyYXX3wxEyZMoG3btqxdu5bvv/8egA0bNtCwYUMeffRRHnvssajOYq677jruvPNOADp27MjYsWMBwpaJHjNmTNhy0K+88goLFizgu+++Y82aNfzjH//gsMMOA2DWrFnMnz+fVq1alf7iKpiktldUdQowpUjcHVHS9kymLNG44grYsAG2bHEv4/jxMGEC/POfUAWen2GkLR9++CEzZ84sXIY6NzeXli1bcuyxx7Jo0SKuvvpqTjjhBI455pi4yotmGgpdJhrCl4P+4osvOOecc8jIyKBZs2YccsghzJo1i9q1a9OtW7cqoQQgjWYWR6N5c3j88WA/P98pg//8Bx54IHVyGUZlpDItTaaqXHDBBdx9993Fjs2ZM4d33nmHUaNGMWnSJEaPHl3m81TGZaMTTVotQx0PN9zgfkePho0bUyuLYRjROeqoo5g4cSJr164F3OiiX375hTVr1qCq9O3bl+HDh/PNN98AbiTQ5s2bEyrDoYceyoQJEygoKGDVqlVMmzYtzFFOVSHtWwRF6dIFevaEqVNhzJhAMRiGUbno0KEDw4YN46ijjqKgoIBatWrx5JNPkpGRwYUXXoiqIiLcf//9AAwcOJBBgwZF7SwO7SMA4lpTqE+fPkyfPp2OHTsiIjz44IPsvvvuCb3OiiBtlqEuDW+9BSedBC1bwtKlNr/ASG9sGeqqR2mXoTbTUAR694Z993XObF59NdXSGIZhJBdTBBGoUQOuv95tP/FEamUxDMNINqYIonDOOZCZCV984XwWGIZhVFdMEUShfv3Aic1rr6VWFsMwjGRiiiAGZ53lfidOTK0chmEYycQUQQxOOMF5NPvqK/jll1RLYxiGkRxMEcSgXj048US3baOHDKNktuzYwrBPhtF0ZFNq3FWDpiObMuyTYWzZsaVc5YoIAwYMKNzPy8ujadOmnOj9QSdPnsx9991XprJPO+00OnXqRLt27WjQoAGdOnWiU6dOfPnll6Uq5+OPP2b69OkRj40ZM4amTZsWlt2pUycWLVpUJnmTgU0oK4Ezz3SmoVdescllhhGLLTu2cPCYg1m6finb8pzDwbV/rmXElyOYtGAS0wdNJ6t2VpnKrlevHnPnziU3N5e6devywQcf0Lx54Ofq5JNP5uSTTy5T2W+88QYAU6dO5YEHHuCtt94qUzkff/wxTZo04eCDD454vH///jwcurplEYouV62qqCo1apT8vZ6fn09GRkbphfawFkEJ9O7tWgYzZ8JPP6VaGsOovIycNjJMCfhsy9vG0vVLGTltZLnK7927N2+/7ZwZjh8/nnPOOafw2PPPP8+VV14JwPnnn8/VV19N9+7d2XvvvXmtHKM9Zs6cyeGHH06XLl04/vjjWbVqFeAWqMvOzqZjx44MGDCApUuXMmbMGEaOHFmq1sSHH35Iz549OfHEE+nQoQNLliwhOzub/v370759e1auXMm4cePo0KED++23H7feeivglEbDhg259tpr6dixIzNmzCjzNYIpghKpWxf8Dw3rNDaM6Dw+6/FiSsBnW942nphVvkk5Z599NhMmTGDbtm3MmTOHgw46KGralStX8sUXX/DWW28xZMiQMp1v+/btXHPNNUyaNInZs2czYMAAbr/9dgBGjBhBTk4Oc+bM4bHHHqNt27YMGjSIwYMHk5OTQ/fu3YuV99JLL4WZhnbs2AG45aoff/zxQp8FCxcu5LrrrmP+/PmoKkOHDuWTTz7h22+/Zdq0aYUtlo0bN3LYYYcxZ84cunXrVqZr9DFFEAf+6KHXX0+tHIZRmVn357rYx3NjHy+Jjh07smzZMsaPH09vf2x3FE499VRq1KhBdnZ24Vd8aVmwYAHz5s3jqKOOolOnTtx3330sX+7csLdv354BAwbw0ksvUSvONWj69+9PTk5OYfDXOiq6XHXbtm0LF677+uuv6dWrF02aNKFWrVr069ePzz77DIDatWtz2mmnlenaimJ9BHFw5JGQkQGzZ8PWrc5UZBhGOI13aczaP9dGP163cbnPcfLJJ3PjjTcydepU1q2Lrljq1KlTuF3W9dRUlY4dO0b0Xfzee+/x6aefMnnyZP75z38yZ86cMp0Dyr7Mdd26dZEEeQqyFkEcZGU5T2b5+ZDk9e4Mo8pyedfLyayZGfFYZs1MLut6WbnPccEFFzBs2DA6dOhQ7rJKIjs7m19//bXQ/r5jxw7mzZtHfn4+K1asoFevXowYMYK1a9fy559/JmWZ64MOOohPPvmEdevWkZeXx4QJEzj88MMTeg4wRRA3vsmvlCPKDCNtGNxjMG0btS2mDDJrZtK2UVsG9xhc7nO0aNGCq6++usz5o7mpjESdOnV47bXXuP766+nYsSMHHHAAX3/9NXl5efTr14+OHTvSuXNnbrzxRurXr88pp5zCxIkTOeCAAyJ2FhftI/j6669LlKFFixbcfffd9OzZk06dOnHwwQdzwgknlOqa48GWoY6Tl192voxPOgkmVwpvy4ZRMZRmGeotO7YwctpInpj1BOty19G4bmMu63oZg3sMLvPQUaP0lHYZausjiBO/U/7LL527vmQ58TaMqkxW7SzuOuIu7jrirlSLYpQCMw3FSevW0KwZrFsHixenWhrDMIzEYYogTkSCVsFXX6VWFsOoaKqaCTmdKcuzKlERiMguInK7iDzt7e8jIieWQb4qj3UYG+lIZmYm69atM2VQBVBV1q1bR2Zm5NFb0Yinj+A5YDbgT137FXgVKNuCHFUYUwRGOtKiRQtWrFjBmjVrUi2KEQeZmZm0aNGiVHniUQRtVfUsETkHQFX/lETNYqhidO7sHNnPmwcbN0KDBqmWyDCST61atWjTpk2qxTCSSDx9BDtEpC6gACLSFtieVKkqKZmZ0KWLGzUUxxBgwzCMKkE8iuBO4F2gpYi8BHwE3JxMoSozZh4yDKO6UaIiUNX3gdOB84HxQFdV/STJclVaQucTlIV16+CJJ2BL+fx0GIZhJIx4Rg19pKrrVPVtVX1LVdeKyEcVIVxlxG8RfPopzJ9furwFBdC3L1x+OQwalHjZDMMwykJURSAimSKyG9BERBqJyG5eaA00j5avurPXXjBwIOzYAeefD3l5kdPt3AlFV799/HH4xGtLvfJKcfeXK1e6he0MwzAqklgtgktww0b/5v364U3gseSLVnl56CFo2dJ5Lbv//uLH58yBDh1gzz3hjjucsliyBG72elb69HG/l18Oq1fD9u1w9dVOyZg7TMMwKpoSF50TkatU9dEKkqdEUrXoXFE++ACOOcYNJ505E/bf340meuopuPZaV7n7HHaYa0FMnw79+sG4cS7vhx+63zVr4NtvXdrmzWH5clvLyDCMxBJr0bm4Vh8Vkf2AbKBwupqqjk2YhKWgsigCgMsugyefdP4KGjRwX/6+OWjQIDjtNLjwQvj9dxe3554wdy7sthv88otrNWza5I7tvTds2AB//OH6HuJc7NEwys3gwe69ff55iMNPulFFiaUI4uksHgY86oUjgBHAyQmVsIoyciTsu68bAfTrr+7PtOuubsnqp592ju+/+w6OP97NQXj2WacEAFq1cn0GtWo5V5jffAPHHeeOffBB6q7JSC82bIAHHoAXXyz94Aej+hDPzOI+wP7At6o6UET2AMYlV6yqQVaW6w9YvTqIa9zYObz32X13mDLFmYpCvOcBzr9Bnz5B/NFHOyXywQeuz6Ak3nwT3nkn2G/TBm66ycxKRvyEToz86ivYb7/UyWKkjngUQa6qFohInojsCqwGWiZZripD7doQz7IeRZVApPijj3a/U6e6PgXPt3VEVq92LYntReZ4H3QQ9OxZsjyGAeHzYb78Ei66KHWyGKkjHovgLBFpCDyNGzX0DWALMSeB5s0hO9uZmqZPj5328cedEujRw01QO9FbD/b995Mvp1F9KKoIjPQknpnFl6vqBlV9EjgaOE9VB8ZTuIgcJyKLRGSJiAyJcPxSEfleRHJE5AsRyS79JVQv/FZBrH6C3FwYNcpt33svXHopXHml2zdFYMRLfn5gGqpVC374AdauTa1MRmoo1RgBVV0GbPN9E8RCRDKAUcDxuBFH50So6F9W1Q6q2gnXCf1gaeSpjsSjCMaOdX/YLl3c0FSAQw91pqRvvnHLWBhGScybB5s3O+97/oz5klqiRvUk1szijiLyvojMFZF7RGRPEZkEfAzEM77gQGCJqv6oqjuACcApoQlUdVPIbj28FU7TmcMPD+YmrF9f/HhBATzoqcsbbww6hnfZBQ45xM1l+ChtFwApmWXLXMU3pFj7tGoxdqzrXxIpHmrVghEjSi7DNwV17x55Da1//Qvatg2fAZ+XB7ff7t43/3x16sBzzyXu2owUoKoRA/A1bqG5fYFrgFW4r/bMaHmK5O8DjAnZ/z/gsQjprgCWAsuBfaKUdTEwC5jVqlUrre4cdpgqqE6aVPzYm2+6Y61aqe7cGX7sX/9yxy68sGLkrIpcdpm7R6D67ruplqZs7Njhnr9/HZFCRobqrFmxyzn3XJf20UdVJ09224cf7o6tWqVap05Q3iWXqP7wg+ohh0Q+X/Pmqtu3J/3SjXIAzNJo9XXUA5BTZP/HaGmj5I9LEYQc7we8UFK5Xbp0SdZ9qjTcfbd7Mv37qxYUBPG5uaoHH+yO/fvfxfPNmhUoCT/f7NlOMfz+e8XIXplZs0Y1MzO88lq/Prnn3LhR9dJLVT/+OHFlvvyyk/+vf1XNy3PPOjRcc4073r696rZt0ctp186lmz1bdfVqt73LLk7RDBsWnKN27fBKf889VT/5xJ0rP181O9vFv/hi5PPk56vefrvqM88k7h4YpaesimAhcADQ2QsLQvej5QvJ3w14L2T/FuCWGOlrABtLKjcdFMF33wV/ur59VTdsUF24UHX//V1co0augilKfr5q48YuzaJFLk3LlsEXXbozfLi7F8ceGyjU885L7jnvuCN4Zr/+Wv7yCgpUu3RxZT71VOQ0W7eq7rOPS3PzzZHTrFoVVPx+y9LP8/nnqk2auO2pU1W/+SY4dvzxTmmE8swz7linTuEfLj4PPRS0Un75pezXbpSPsiqCT2KEj6PlC8lfE/gRaAPUBr4D2hdJs0/I9kmxBPVDOigCVdXx41WzstwT+stfVOvVc9vt2rkvuGicdZZL99hjriXgK5TMzOJ/4NLy3XeuUijKypWq771XvrITTV6e6ttvBwozN1d1993dvfjoI6dY/dbBm2/GX+7y5a5yjIc//wwqVFA94YTIFWVpmDrVldW0qSs/Gl9+qVqjhgvTphU/7psYe/YM4s47L6jQwSkcX94tW1S/+sp9bBRl2zbVPfZweT78MPxY6H0G1RtvLPUlx7KbIDQAACAASURBVEVBgbump56KHCK9t9FYs0b1ueeCvC+8EPteF+Wrr1Tnzy/1JSSdMimCRASgN/CD1wdwmxc3HDjZ234EmAfkeAqmfUllposiUFVdvFi1c+fgT9SvX+SWQChjxri0rVu73zp1VLt2ddt33VU2OfLy3JetSHHb89atznzgV7CVhdtuczLtvbfqjBmqTz/t9g84IKjc/C/Vdu3iq6DHj1etXz9+5fHkk1poomnY0G0/+2z5ruvEE105w4aVnPamm1za+vWd7KHcfLM7duutxeX1Q9E8sfDNmccfH8Tl5QUtL79vYdddS36Hy4L/3kcLdeu6Po6S+PBD1WbNiuc/99z45Pj8c/c/2XVX1Z9/Lt81JZqUKYJkhHRSBKrua2vECNUJE+KrrH7+OfwFHjHC2XPj+Yr0yc1V/e03FxYscB2IoWWG2p59ezSoXnBBea40cWzeHFS8oFqrVtAaGDcuSJeXF/zpi3as7twZ3IPly1Uvvjj8HvToEVuG/PxAQY4frzp2bFARzpkTlB0aduyIXeaCBYFyX7Wq5PuwbZszLfoyDxrkTDO//abavbuLe+utIP2cOUHaSIMRYrF2ratsQfWLL9w5fOXg98X4gyAefDD+cuNh2bJAQfftq3rRReHBV0bdu7tn7lNQ4Fqzv/3mzHa33+4qcVA96CCXd9Cg6C3H3Nzw/S1bVNu2De7h0UeXvwWYSEwRpBl+BeS/+AUFQcti9OjYeceODf5UoaFZM1dphNqefTNFjRpaaAevDCNHHnnEydOtm+pVVwXX0KJF8cr2yivdsZtuCuLy8lzLoeg9qFPHVWINGrj96dOjy1B0dFdBgeoppxQvMzS0bev6gyKxcqXqoYe6dBddFP+9KChQfeKJ8BFAoWHt2vDr3nVXFx9pMEJJhI7ICg3vvBP5niSC/HzVI4905Z52WuSK948/XAc3qI4c6eIWLAhMYKFBxLW2QhWG33LcYw9nNsrPdx9YtWur9uoV9P3471L79kFf3RNPJOY6E4EpgjTjuefcV9DixUHcSy+5p73vvpHtvJs3BzZiv/XQrJkLZ5wRfIGG2p79P9cdd6jut5/bfvvtirjC6OzcqdqmjZPljTdc3BtvqP7jH6qvv148/WefaaEpza9EJk0KKn7/HnTrppqT4477ZpW+faPLEenr9/ffnZnOLzM0+H1AAwcWL+v994MWTdOmqj/+WPr7kpPjPgxCzxlJodx7r+tM37Sp9Of4+Wf3Hvjl77mnM9H5FG0lJYJRo1x5TZrEbiW9/XbwTO+9N7jfWVmBvPvvH3l0V35+8DxPOcWZv0KVR5MmwUCEmjVVv/1WdeJEt1+vnurSpYm51vJSbkWAc03ZHTjMD/HkS0YwRVA2duxwX8Tg7OWhfPed6t/+poW21Kefjt2k9W3P4L6qtm8P/giho3CmTXPNY78CTRSjR7uOzM6dXTjmmKAD3f8DtmsX/lUXjfx81b32cnm+/trF+WaTRx6JnGfFCveHr1EjcqX81ltaaAaK1x4+f37w1f6//7m4nTudDd83V/Tq5cwYVRm/H6JDB9V16+LL89JLThH7zzs0+GabV18tuZwLLgivwOPpc/NZujRQHqC6226u9XzUUeFlDh8e5PEHbuy5ZyBvv34lmwCTRbkUAXA/sAyYAvzPC5NLypesYIqg7PznP8ELe8UVzsb55JNBBZSdrfr99yWXk5vrvp522cUpEVU3OgSc2WTbNmcTbt5cC1shpRl1EYvt28Pt/36oXdtd30EHuf1Ro+Iv0+/nuOEG1+IBd47Nm6Pn+b//c+muvjqI27lTdejQoOIO/RqOh5EjtdAMl5Pj+iF809vw4fEptsrOn38G70WrVpFHNPls2eJaSLHMaeDm28TDxo3uA6FuXTfktbT2e3/AwSGHBMNg8/NdCyMjwymr0Ep+zZrgIyM0RJtvkWzKqwgWAXVKSldRwRRB2SkoUH34Ydd56psZ/JfzwgvdCKB42bKleFPcn+cweXK4mQlUr78+PO/o0ZFHVbz+evBlHgn/a7t9e9fBO3Om6uWXh5+rcePSXcu0aS5fy5aqp5/utocMiZ0nJ0cLm/633+7MY37FLeL2S2sHz8sLyvDDXnvFP1y1qvDjj6oHHuiuLyPDPb877ige/v53lyYz0w2HnjWreMjJKd193rKlfKOWVq6MbFpdsyZy/9gffwSy/vOfWtiKTkUncnkVwTtAVknpKiqYIig/M2e6YZW+jfSllxJTrv+i77tv8Ad+8UX3Zxdx9vg5c4I/+F//Gl5hv/ZaUJFHaz77yyLcc094/GuvBZ24Q4eWTu78/GDiHThFGc/kr6OPDq+0fTNAeWYR//BDMPqmd+/yz/2orGzf7lpgJX3t//3v8bVSqwKhc1mKzreoCGIpgnic10/CeSj7CCh0g6KqcfjQSjyVyWdxVWbTJnjpJTjmGLewWCJYuhTatQv2H3wQrrsOhg51y2XvuadbSG/btiDNtdfCQw85Rzvt2wfLIL/zTuC602fbNthjDyf7Dz/APvuEH//5Z/jkE+jXL7ZTn0jccEOwmN+558ILL5ScZ/ly5+Jx5063X6+ey7v77qU7d1FmzHA+rU8/vfr7EP70U+eIKRKNG8PAge6+VhfuvhvuuMO5r50yxcXl5cHKldCyiLsvVVi0yPkn8WnUqOz/11g+i+NpEZwXKZSUL1nBWgSVG3/5g0MPDZrQ27erduwYfOVdcIFrHfgthU8/DUwyfuff+ecXL9sffnjAAYmX++uvA/n8fg/DSDRr1gQtvrlzXQvQH9p96aVBX9r69W60XtEW0llnlf3cJGDUUG1gPy/UiidPsoIpgsrN+++7IXbLloXHL1jg4l9+OYjz1+HxTTr167sVQf24ojbX/v3dsX/9K/FyFxSoXntt+KgPw0gG/nyLgw4KlpHxQ8eObvKovzJAVpb7uPLDLbeU/bzlUgRAT+Bn4FPgM+AnbPiokQC2bw+f1OMPa+3Qwe37wyhV3ZeS/6epLOOyDaMs/PBDMLIMVM880w0I8Cdr+qFLF9UlSxJ33liKIB4L5L+BY1T1cFU9DDgWeKgsNirDCKV2bWeLb9AATjsNLrzQxZ91lvudODFI++67zlbatSvsvXfFy2oYiWKffeDii6F+fXjqKZgwwTmkmj0b/u//ICPD9Z1Nm5a4/ruSiKezeI6qdiwprqKwzuLqx/btzquW3zGaM38LB1wxEg58HNllHY13aUzjHy9n0bODGXFPFoMHp1Zewygvqs7bYEZG8WPbtzuvb4kmVmdxzTjyzxKRMcA4b78/zluYYSSE0Jd+y44tDJh6MHLoUjRjGwqs/XMtaxuPgEGTOOG06UBWqkQ1jIQgElkJQHKUQEnEYxq6DOej+GovzPfiDCPhjJw2kqXrnRIIo9Y2au2xlFeWj0yNYIZRjSnRNFTZMNNQHBQUuAH7jRunWpJS03RkU9b+uTb68V2asnrw6vgK27ABGjZMkGSGUbWJZRqK2iIQkYne7/ciMqdoSJawRjkpKIAzz4RmzWDmzFRLU2rW/bku9vHc2McLGTbMzb4ZNMjN2DEMIyqxTEPXeL8n4txIFg1GIpk9202bveee+PMsWOCmZYa26u6/HyZNcpXfY48lXMxCFi6EDz5IeLGNd4ndimlcN45WzmOPwfDhbvuZZ6Bv3/DpzIZhhBFVEajqSm/zclX9OTQAl1eMeGlCfj5ccolbZ+GOO9y4sZLYsgUOOQSOOMKNOdu8GT7+2K3n4PPqq7BxY+Llzc2FXr3c+hRvvBF+bPt2mDXLtUxCyc+H6dNd3hhc3vVyMmtmRjyWWTOTy7qW0D316qtwtbf6yc03O9PQf/8LJ5zg7pFhGMWIp7P46AhxxydakCrLjh2u8isPTz/tWgQ1ariv+/PPh61bY+d59ln44w+3/dJL0KULnH22q4CHDnUDk3Nz4ZVXyifbtm3Fr+/pp93iKACXXhosEJSbC0cdBf/4Bxx0kFtIRtUtHNSpE3TrBvvuC+PGFVcUHoN7DKZto7bFlEFmzUzaNmrL4O43wujR0L+/++qfPNmZwJ5/Hq65BgYMcOf85z/hvvucDM2aOSXpT1QwDCOcaDPNcCODvge2AnNCwk/AuGj5kh0q1czinTvdNNi2bcu+4P6aNc7HIziHuv602quuin1e3w3X/fcHecB5ysjLC5zkHnhg2eTautU5nc3Kco4K1q938bm5wSLrvqebvn2dTKeeGj410l9i1N/2FxICt8BKlGU6N2/frHd8fIc2HdFUa9xVQ5uOaKp3fHyHbv5tmepJJxU/R9Fw1VXh6/z+8INzngCqH3xQtvthGFUcyrLEBNAAaA2MB/4SEnaLlqciQqVSBK+/HlQ+oYvoqDrPFePGuYXIc3NdmDhR9YQTnG+7M85wC/f7njd8T9fffOPcX4GriEeNcmHBgqDsom64/vxT9brrVI87Lli3eOvWwAFtpHV8d+50PhxXrAiPz8tTff75wHuIH045xcn3+ONauCjK0qXBug++W6+GDd0KbnfdFbh0athQ9YEHnEzPPRfureOkk5x7rljk56tOmRIonoYNndPYG290im+//dw8/Xvvdev7Rlrs/V//cnn33bd8jpVXrHBOESqTV3LDiIMyKYJiCWF3oJUf4s2X6JAyRbB0aXGPFqF+6o48MojPzw//Ss/ICPdzVzTUquVcfPnceWfxNJmZbvGdgoL43XBdeqlLd+21xY9dcYU7VreuW8B/0yb3tRy6+E/nzs5Zge8S7O67g4X7fd+Avu9BX8YvvgjO8dtvruIP9ZCuGt7a8O/PUUepDh7slOeHH7rw/vvOV2Oos4CDD1b96aeSnlZxtm8PHCXcd5+LKyhwSshv7ZTEunXBamATJ5ZeBsNIIeVSBLgRQos9E9FPQAEwr6R8yQopUQSLFztfiO3bByYg3zdj3bqBycN3YOt7P2/UyHnWqFHD7R9wgHMR9v33rjLyPbTcdVf4+XbscBXlZZe50Lt3UGFed53b3m23kt1wzZzp0jZu7PxH+owbF5TnV7ChiqplS+dRxl9HevLkcKWUnR0cKyhwJqHMTNX//rd093XlStVLLgnuT6zQurXzfFMeh6/vv+/K2mUXp2x9R83t25fcSsjPd605X54OHSK7qjJSx86d1cOfZ5IoryL4DmgMfOvtHwE8U1K+ZIWUKIK77goqgBtucHG+o9tBg1QHDHDbt9/uKkb/q/qxx1zaP/+M7PKqoED1999LNjMUFLj1Z0Mrxngc4hYUBI4ATj/dKarvvw/s5U884fw0Hnyw269f35lQIvV3hJ6/qBmsoKB8/v9+/dUpkTvvdHL26hWEQYPc0oyJqnT79o2saPxWgs///uee3++/u/177w2U+557uu3XX0+MTEb5WLPGOY+uVcuZHW+6yS32b4RRXkUwSwOFUMPfLilfskJKFEF2dlBhiKi+916wiP6336p+8onbbtHC2d3BVRa5uYmVw/dunpnpvqbjYfLkwDt97drOMzo47+u+AioocAphzZro5ezc6RwR9+1btb+6Vqxw3sf79HGV/ZQpQSvB90g+YUKwTrBvtqpRw8VNmaL66KPuWDzOZ195xbVmLrmkdI6UjZLZudP9J/z/YtFwzjnWaguhvIrgQ9wqX496HcePAF+WlC9ZocIVwdy5WmiKGTw4qFBBtVs3lyY/P3ACvNtu7vfhh5Mjz7RprgO6NCxbFrRawHWuWqUU0KePuy99+qh+9JH7svT7I/yOe3CedFSdgvdbBZMnRy7zzz9d5R9aMbVvb1+qiSI3N3yU2tFHO0/2X3yhetFFgamzqHPrykBBgevnWry4Qk9bXkVQDzffoCbOTeXVQOOS8iUrVLgiuP12LTQBbdsW3gk8blyQ7p57gvg99qicFe2sWapDhqj+/HOqJalc/PJLYC7z/Qhee637w65apfqf/7hRSqEtoYcecum6dg1vFeTmulaA/57Uru3cnvn9EXXruhZHvBQUlG+UU0WyerXqlVe6e7NqVenyLl/uzHODBwd9baru+t96y5U7frxTsBs3qh5xhLufDRu640WZMsW14ERc31Cqyc93gzsOOyxowYhE/5BIAuVVBG2AzJD9ukDrkvIlKyRVEeTmuj+9X1EWFAQjTfyX6ZtvAltkaAfsihVBp+cDDyRPRiM53HdfoMjPPrtkk8LWraq77+7SN2umeuyxrtUVaqb461+d6VBVdcsW54jZVw5ff12yTOvXu36SrCzVZ54JVzjz5rl31S8/NH7ECNVnn3Xvaug7mkx+/jn4r4BrSZ10knNOHUpBgetbue02F2691X3Nh7rsql1b9frrXd9Qr17hraoGDdy8Hf++x3Iw7ftCbdIkMPulgjVrggEffvBH4jVoUGEu98rdRwDUDtmvDcwsKV+yQlIVgW/7/ctfXAdmTk7wIu3cGaSbPz/yV/Utt6gec4z70xtVi+3bnWlo4MD4K89JkyLbpzt3dqbBzZuL57nySpemVavYfTIrV6ruv394uf36uclxRUdadeyoevPNzrdhUVlq1nR9HC++mLz3csGCYI5Hhw5OAfgj0kRchZyX54Yo9+tXXEa/8u/TJ3BMHRoaNXKDNLp2DeL23rvkCjQvz/0fwSmPf/876Pwvibw81xfYv79T6DffrLphQ/z3pKDAfRxOnBjMyWnUyE30XLkyGG3n9zWVdUJqKSivIsiJEFc9O4tDhwd27BiMtb/kkuSd06ja5Oe7Cun114OhwbHYvj2YB3LMMU7pvPmmmxB37LFuxMtzzwVfvX/9q+ojjxSfh5KRoXriiUGfVOgX8/nnq551lvtCD/3SzspyrYXysmGD6qefuhbJhRcGMvToEczJWLXKfRj55z/0UDcBEty1DBnihkjffbdzVr1uXVD+7NmuJVCnjmsZhB6bO9d9sMVboa9ZE+4MOCPDVe7RBnLMm+eeQeikRz80aeLOHWkI89Klqi+84OTt1csN2Q7N27178Y/H9euD5zxwYHiLLz/fmZ6PPtq1NEeOdPN8/AmjZaC8iuAD4OSQ/VOAj0rKl6yQNEWwfXvwZ/OXb/BDlKUQDKNM/PKLq1T80UqRvpD9loVva1+4MBgKfPLJwUzzbducErrmGtf3UPTL8o8/3KS/bt2Ccj/8MLps+fmqixa5fo5bbnFf8H446aRgQl3R0Lt35H6xjz4KRqr5H1ihs+RjEdoKLw/btqm+9pq7b37n/xlnhPf5vPlm8RZV27Zu6PiUKW6kWWjHdGhr77nnwufkhLZkevZ0y8BEm/+SkxPMQ/LP99//Rm7d+R8PZaS8iqAtMB34BVgOfAm0KylfskLSFMHUqe52ZGe7Hn1/VMgee1Tt4ZJG5eT994Ov5exsV1lMmqQ6bJhbzuOyy4rPzdixo2yzqn2GD9dCs1TRshcudBW/b+KJFerUcWaaCy90rYLPP4/dp7JqlRt6PHhwhZhAYvLdd4E57/LLXcvAN9eBW5bloovc6KPQL/SCAjc03O8XOvhg11L597+DvCee6Fo4kye7r/94lyF5/fXILZC99nJzfZ56ysnao4eba1NGErXERBaQFW/6ZIWkKYLbbnO345pr3H5Ojvt6efTR5JzPMGbPdh26FbVu0Y4dwZfmoEEubsaMYASOH5o1c5Xabbe5dafGjXNh4kRnOknUl3qqmDo1mFvjV8C1arlBHiUpqh9+cH2IECgFcOa78pCXp/ruu67l9fe/O6Wd4D6dWIogqqtKERmgquNE5PpIx1X1wYgZk0zSXFUedBDMmAFvvw29eye+fMOoDMybB507u+XTjz0W3nvPxdev7xz4nHee83NRI54V6qswkya561WFvfd2y7V3jejFsTi//up8ccyf7zzQP/ec8wlSyYnlqrJmjHy7eL/1Ey9SJWP9eudMpVYtt46/YVRX2reHu+92Tnveew/q1HF+HG65Jb38O59xBrz+uvNlcdNN0KBB/HmbN4fPPnP+Lo47Do48MnlyVhCxFEFb73e+qr5aEcKkjI8/do5SDj0U6tVLtTSGkVxuuAGWLYOdO+G226B161RLlBpOPdWFstC4MYwcmVh5UkgsRdBbRIYAtwDVWxG8/777PTqSMzbDqGZkZMDjj6daCqMSEcsQ+C6wHugoIptCwmYR2RRP4SJynIgsEpElnlIpevx6EZkvInNE5CMR+UsZr6N8+E7YjzkmJac3DMNIJbGc1w9W1YbA26q6a0ior6q7llSwiGQAo3D+jbOBc0Qku0iyb4GuqtoReA0YUeYrKStLl8JPP0GjRq4TzTAMI82IZRoCQFVPKWPZBwJLVPVHABGZgJuMNj+k7E9C0k8HBpTxXKVj7Vp46in45hvXWQSuwycjo0JObxiGUZmIqghE5AtVPURENgMKSMhhjaNV0Bw3Ac1nBXBQjPQXAu9EkeVi4GKAVq1alXDaOLj3Xnj44WC/bl244ILyl2sYhlEFiaoIVPUQ7zfpw0dFZADQFYg4dlNVRwOjwc0jKPcJ1651v1ddBRddBH/7mxs6ahiGkYaUOGtERNqKSB1vu6eIXC0i8Qw4/hVoGbLfwosrWv5RwG249Yy2xyd2OcnNdb+HHgodOpgSMAwjrYln+uAkIF9E2uG+ylsCL8eRbyawj4i0EZHawNnA5NAEInIA8BROCawuleTlwVcEdetW2CkNwzAqK/EoggJVzQNOAx5V1cHAniVl8vJcCbwHLAAmquo8ERkuIid7yUbi1jB6VURyRGRylOISy7Zt7jczs0JOZxiGUZkpcdQQsFNEzsG5qTzJi4vLlqKqU4ApReLuCNk+Kk45E4u1CAzDMAqJp0UwEOgG3KuqP4lIG+DF5IqVZEwRGIZhFBLPPIL5OIf1iEgjoL6q3p9swZKKrwjMNGQYhhHXqKGpIrKriOwGfAM8LSIpWYI6Yfh9BNYiMAzDiMs01EBVNwGnA2NV9SAgNbb9RGGmIcMwjELiUQQ1RWRP4EzgrSTLUzGYacgwDKOQeBTBcNwQ0CWqOlNE9gYWJ1esJGOmIcMwjELi6Sx+lRB/BN4icmckU6ikkp/vHHKIQO3aqZbGMAwj5ZSoCEQkE7cgXHug0JaiqlVzlbbQ/gGR2GkNwzDSgHhMQy8CzYBjgU9xawZtTqZQScVmFRuGYYQRjyJop6q3A1tV9QXgBGIvJ125sRFDhmEYYcSjCHZ6vxtEZD+gAbB78kRKMqYIDMMwwohnraHR3ozi23Grh2YBd8TOUomxoaOGYRhhxDNqaIy3+Smwd3LFqQBs6KhhGEYYsVxVXh8ro6pWzWUmzDRkGIYRRqwWQdJdVKYEUwSGYRhhxPJZfFdFClJh2PBRwzCMMKKOGhKRkSJySYT4S0TkvuSKlUSsRWAYhhFGrOGjvXA+iovyNHBicsSpAEwRGIZhhBFLEdRRVS0aqaoFQNVdm8FMQ4ZhGGHEUgS5IrJP0UgvLjd5IiUZaxEYhmGEEWvU0B3AOyJyDzDbi+sK3AJcm2zBkoYpAsMwjDBijRp6R0ROBQYDV3nRc4EzVPX7ihAuKZhpyDAMI4yYM4tVdS5wnr8vIs1U9fekS5VMrEVgGIYRRjyLzoUyJSlSVCSmCAzDMMIorSKouqOFfEwRGIZhhFFaRfB0UqSoSKyPwDAMI4wSFYGItBWROt7ufBG5WkQaJlmu5GEtAsMwjDDiaRFMAvJFpB3wFNASeDmpUiUTUwSGYRhhxKMIClQ1DzgNeExVBwN7JlesJGKmIcMwjDDiclUpIufghpG+5cXVSp5IScZaBIZhGGHEowgGAt2Ae1X1JxFpA7yYXLGSiCkCwzCMMOJxVTkfuBrA811cX1XvT7ZgScNMQ4ZhGGHEM2poqojsKiK7Ad8AT4tI1XRTCdYiMAzDKEI8pqEGqroJOB0Yq6oHAUclV6wkYorAMAwjjHgUQU0R2RM4k6CzuOrim4ZMERiGYQDxKYLhwHvAElWdKSJ7A4uTK1aSUA1aBNZHYBiGAcTXWfwq8GrI/o/AGckUKmls3+5+a9eGGqVdXcMwDKN6ElURiMhNqjpCRB4FIrmsvLqkwkXkOOARIAMYo6r3FTl+GPAw0BE4W1VfK6X8pcP6BwzDMIoRq0WwwPudVZaCRSQDGAUcDawAZorIZG84qs8vwPnAjWU5R6mxoaOGYRjFiOWh7H/e7wtlLPtAXL/CjwAiMgE4BShUBKq6zDtWUMZzlA5rERiGYRQjlmlocqyMqnpyCWU3B5aH7K8ADopftDBZLgYuBmjVqlVZinCYIjAMwyhGLNNQN1xFPh74mhQ6pVHV0cBogK5duxbrr4gbGzpqGIZRjFiKoBnOvn8O0A94GxivqvPiLPtX3JLVPi28uNRhQ0cNwzCKEXUMparmq+q7qnoecDCwBJgqIlfGWfZMYB8RaSMitYGzgZjmpqRjpiHDMIxixBxMLyJ1ROR0YBxwBfAf4I14CvZ8GFyJm4y2AJioqvNEZLiInOyV/w8RWQH0BZ4SkXhbG2XDTEOGYRjFiNVZPBbYD5gC3KWqc0tbuKpO8fKHxt0Rsj0TZzKqGMw0ZBiGUYxYfQQDgK3ANcDVIoV9xQKoqu6aZNkSj5mGDMMwihFrHkH1W4PBTEOGYRjFqH6VfSzMNGQYhlGM9FQE1iIwDMMoxBSBYRhGmpNeisD6CAzDMIqRXorA+ggMwzCKkZ6KwFoEhmEYhaSXIjDTkGEYRjHSSxGYacgwDKMY6akIrEVgGIZRSHopAjMNGYZhFCO9FIG1CAzDMIqRnorA+ggMwzAKSS9FYKYhwzCMYqSXIjDTkGEYRjHSUxGYacgwDKOQ9FQE1iIwDMMoJL0UgfURGIZhFCN9FEFengs1akDNWB46DcMw0ov0UQShZqHA/7JhGEbakz6KwMxChmEYEUkfRWAdxYZhGBFJP0VgQ0cNwzDCSB9FYKYhwzCMiKSPIjDTkGEYRkTSTxGYacgwDCOM9FME1iIwDMMII30UgfURGIZhRCR9FIGZd0dY6wAABkVJREFUhgzDMCKSforAWgSGYRhhpI8iMNOQYRhGRNJHEViLwDAMIyLppwisj8AwDCOM9FEEZhoyDMOISPooAjMNGYZhRCT9FIGZhgzDMMJIqiIQkeNEZJGILBGRIRGO1xGRV7zjX4tI66QJY6YhwzCMiCRNEYhIBjAKOB7IBs4RkewiyS4E1qtqO+Ah4P5kyWOmIcMwjMgks0VwILBEVX9U1R3ABOCUImlOAV7wtl8DjhRJkh9JUwSGYRgRSaYiaA4sD9lf4cVFTKOqecBGoHHRgkTkYhGZJSKz1qxZUzZpsrKgcWOoV69s+Q3DMKopNVMtQDyo6mhgNEDXrl21TIW8/HIiRTIMw6g2JLNF8CvQMmS/hRcXMY2I1AQaAOuSKJNhGIZRhGQqgpnAPiLSRkRqA2cDk4ukmQyc5233AT5W1bJ98RuGYRhlImmmIVXNE5ErgfeADOBZVZ0nIsOBWao6GXgGeFFElgB/4JSFYRiGUYEktY9AVacAU4rE3RGyvQ3om0wZDMMwjNikz8xiwzAMIyKmCAzDMNIcUwSGYRhpjikCwzCMNEeq2mhNEVkD/FyKLE2AtUkSpzJj151+pOu123XHx19UtWmkA1VOEZQWEZmlql1TLUdFY9edfqTrtdt1lx8zDRmGYaQ5pggMwzDSnHRQBKNTLUCKsOtOP9L12u26y0m17yMwDMMwYpMOLQLDMAwjBqYIDMMw0pxqrQhE5DgRWSQiS0RkSKrlSRYi0lJEPhGR+SIyT0Su8eJ3E5EPRGSx99so1bImAxHJEJFvReQtb7+NiHztPfdXvGXQqxUi0lBEXhORhSKyQES6pcPzFpHrvHd8roiMF5HM6vi8ReRZEVktInND4iI+X3H8x7v+OSLSubTnq7aKQEQygFHA8UA2cI6IZKdWqqSRB9ygqtnAwcAV3rUOAT5S1X2Aj7z96sg1wIKQ/fuBh1S1HbAeuDAlUiWXR4B3VfVvwP6466/Wz1tEmgNXA11VdT/c8vZnUz2f9/PAcUXioj3f44F9vHAx8ERpT1ZtFQFwILBEVX9U1R3ABOCUFMuUFFR1pap+421vxlUKzXHX+4KX7AXg1NRImDxEpAVwAjDG2xegF/Cal6TaXbeINAAOw/nzQFV3qOoG0uB545bOr+t5NNwFWEk1fN6q+hnOR0so0Z7vKcBYdUwHGorInqU5X3VWBM2B5SH7K7y4ao2ItAYOAL4G9lDVld6h34E9UiRWMnkYuAko8PYbAxtUNc/br47PvQ2wBnjOM4mNEZF6VPPnraq/Ag8Av+AUwEZgNtX/eftEe77lruuqsyJIO0QkC5gEXKuqm0KPeS5Aq9VYYRE5EVitqrNTLUsFUxPoDDyhqgcAWyliBqqmz7sR7uu3DbAXUI/i5pO0INHPtzorgl+BliH7Lby4aomI1MIpgZdU9XUvepXfRPR+V6dKviTRAzhZRJbhTH+9cLbzhp7pAKrnc18BrFDVr73913CKobo/76OAn1R1jaruBF7HvQPV/Xn7RHu+5a7rqrMimAns440oqI3rVJqcYpmSgmcXfwZYoKoPhhyaDJznbZ8HvFnRsiUTVb1FVVuoamvc8/1YVfsDnwB9vGTV8bp/B5aLyL5e1JHAfKr588aZhA4WkV28d96/7mr9vEOI9nwnA+d6o4cOBjaGmJDiQ1WrbQB6Az8AS4HbUi1PEq/zEFwzcQ6Q44XeOHv5R8Bi4ENgt1TLmsR70BN4y9veG5gBLAFeBeqkWr4kXG8nYJb3zP8LNEqH5w3cBSwE5gIvAnWq4/MGxuP6QXbiWoAXRnu+gOBGSC4FvseNqirV+WyJCcMwjDSnOpuGDMMwjDgwRWAYhpHmmCIwDMNIc0wRGIZhpDmmCAzDMNIcUwSGUU5EpHXoKpGGUdUwRWAYhpHmmCIwjAQiInt7C8H9I9WyGEa81Cw5iWEY8eAt+TABOF9Vv0u1PIYRL6YIDCMxNMWt/XK6qs5PtTCGURrMNGQYiWEjblG0Q1ItiGGUFmsRGEZi2AGcBrwnIltU9eVUC2QY8WKKwDAShKpu9ZzlfOApg2q57LlR/bDVRw3DMNIc6yMwDMNIc0wRGIZhpDmmCAzDMNIcUwSGYRhpjikCwzCMNMcUgWEYRppjisAwDCPN+X/3rrGutsHhdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.datasets import make_classification \n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X, y = make_classification(n_samples=1000, n_features=30,\n",
    "                           n_informative=3, n_redundant=0,\n",
    "                           n_classes=2, n_clusters_per_class=1,\n",
    "                           class_sep=.6, random_state=42)\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)\n",
    "\n",
    "\n",
    "k_values = np.arange(1, 100, 1)\n",
    "train_error = []\n",
    "test_error = []\n",
    "\n",
    "for k in k_values: \n",
    "    knn = KNeighborsClassifier(n_neighbors=k)\n",
    "    knn.fit(X_train, y_train)\n",
    "    # Note: miss-classification error = 1 - accuracy\n",
    "    train_error.append(1-knn.score(X_train, y_train))\n",
    "    test_error.append(1-knn.score(X_test, y_test))\n",
    "    \n",
    "\n",
    "plt.plot(figsize=(10,5))\n",
    "plt.plot(k_values, train_error, \"r-\", linewidth=2, label=\"Training Error\")\n",
    "plt.plot(k_values, test_error, \"b-\", linewidth=2, label=\"Test Error\")\n",
    "plt.plot(k_values[np.argmin(test_error)], np.min(test_error), \"og\", markersize=8, label=\"Min. Test Error\") \n",
    "plt.xlabel(\"k\")\n",
    "plt.ylabel(\"Miss-Classification Rate\")\n",
    "plt.title(\"knn: Train and Test Error with varying k's\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum test error (best hyperparameter value) is achieved when $k=36$:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_k = k_values[np.argmin(test_error)]\n",
    "best_k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Model Selection and Validation Strategies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Holdout Method\n",
    "In absence of a separate test set on which to evaluate the generalization ability of a regressor/classifier, it is necessary to create one 'manually'. The solution comes by (randomly) splitting the data into a train and a test set. This can be done automatically with the `train_test_split` function of `sklearn` (when performing this operation, it is important to set the `random_state` parameter in order to obtain reproducible results). Thus, the training data can be used for model training, and the test data for model tuning/evaluation. This method is known as **holdout method**. \n",
    "\n",
    "However, using only a test set to tune the model and evaluate it can be misleading. For the same reason for which we cannot use the training set to evaluate a model, we cannot rely on a unique test set to perform both tuning and performance evaluation. This is because the selected hyperparameter might be tailored to the specific data points present in that test set, leading to an overly optimistic perfomance measure. In this way, when implemented on new data, the algorithm performance might deteriorate.\n",
    "\n",
    "There is a workaround for this: instead of splitting the data into train/test sets, we can split it into three subsets: the train/validation/test sets. These three sets can be used as follows: \n",
    "\n",
    "1. train the algorithm (with different hyperparameter values in case of tuning) on the training set \n",
    "1. evaluate all the trained algorithms on the validation set; in this way, the validation set can be used to select the best model hyperparameter (and obtain a first estimate of the model performance)\n",
    "1. once the best model has been picked, assess the prediction performance by retraining it on the training+validation set, and testing it on the test set \n",
    "\n",
    "\n",
    "<img src=\"./img/model_selection/holdout.png\" width=\"700\" height=\"100\"/>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Typical split ratios are: between 50-70% for the training set, and between 15-25% for validation and test set. This, of course, depends on the dataset size. \n",
    "\n",
    "While the holdout is recommened in regression settings, it can lead to unbalanced (or mis-represented) classes in classification. In the figure above, we can see that -due to the randomness of the split- the positive (blue) class is over-represented in the training set, and not present at all in the test set (where there are only instances from the red class). This, in turn, will lead to misleading performance metrics. \n",
    "\n",
    "A solution can be found by performing stratified random sampling, where the stratum variable is given by the class. In this way, it is ensured that each class retains its original proportions within each subset. \n",
    "\n",
    "<img src=\"./img/model_selection/stratified_holdout.png\" width=\"700\" height=\"100\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to see now how to perform stratified holdout with the [iris dataset](https://en.wikipedia.org/wiki/Iris_flower_data_set) provided by `sklearn`. Such dataset has 3 classes. Notice, in both cases, the use of the `stratify` parameter. (Normal holdout, without stratification, can be implemented simply by omitting the `stratify` argument). We will perform a 60/20/20 split (note that for the validation set size, 0.2/0.8=0.25). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(150, 90, 30, 30)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.datasets import load_iris\n",
    "X, y = load_iris(return_X_y=True)\n",
    "X_trainvalid, X_test, y_trainvalid, y_test = train_test_split(X, y, \n",
    "                                                              test_size=0.2, \n",
    "                                                              stratify=y,\n",
    "                                                              random_state=1)\n",
    "# validation_size = 0.2 / 0.8  \n",
    "X_train, X_valid, y_train, y_valid = train_test_split(X_trainvalid, y_trainvalid, \n",
    "                                                      test_size=0.25, \n",
    "                                                      stratify=y_trainvalid,\n",
    "                                                      random_state=1)\n",
    "\n",
    "original_size = X.shape[0]\n",
    "train_size = X_train.shape[0]\n",
    "validation_size = X_valid.shape[0]\n",
    "test_size = X_test.shape[0]\n",
    "(original_size, train_size, validation_size, test_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check the class distribution within each split:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAHiCAYAAADbK6SdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3debglVX3v//dHBkVBiNJGZRBUUHHWFufhRrwBBzBXo+DIjUpMJEbl5gYTB8Ro1Bg1RpJIjBfjBKjRtIohP0fUiNI4oICYFgcaRBoEBBQB+f7+WOvg5nCGfbr36drd/X49z3nO3lVrr1q1qr71rWnXTlUhSZKGc7OhGyBJ0pbOZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLApj4ZJ/mLJO+adNkx6qokd51EXZubaeubJI9Mcs4E6/tUkuf214cm+dIE635mkv+cVH1boiRbJbkyye6TLLupS/K+JEcN3Y5RSc5J8sgJ1fXcJJ/qr7fu26E9JlT34OvJRk3GfcP27SS/SHJhkn9MstNCn6mq11fV88epfyllN0SSzye5OskVSX6e5PQkRya5+RLq2CgJbX2mk+QOSf4lyU/6PH43yWuS3Gq52rlAW45Kcm1vxxVJvpfkHUnuMFOmqr5YVXcbs673LVauqg6oqvdMoO179P7feqTu91fV/9zQujclfSM383d9kl+OvH/mUuurql9X1fZV9eNJll2qJL+V5Li+Lft5Tzx/NuZnl5w407w0yZlJrkqyNsmJSe61XjOwAZLcta/bM8vxwiQfT/LY0XJVdbeq+uI4dS02zap6T1UdsKFt79P8UpJDR+petvVkXBstGSc5Angj8GfAjsBDgDsB/1+Sbef5zNZzDZ8Sh1fVDsAdgCOAg4GTkmTYZm2YJLcBvgJsBzy0z+PjgJ2AuwzUrBN6O24D/B5we+D00YQ8CX1jN/VnizY1fSO3fVVtD/wYeNLIsPfPLj/lcT/q7cC2wN1p8fFk4PvLOL1jgD8GXgT8FrA38HHg8cs4zQWNLNf7A58FViV51qSnswmtE+uvqpb9D7g1cCXwtFnDtwfWAX/Q3x8FfBh4H/Bz4Pl92PtGPvMc4EfAJcArgR8C+418/n399R5AAc+lbQAuBv5ypJ59aUnnMuAnwDuAbUfGF3DXeebn88DzZw3bHfgF8MTF6gdO6fVf1fvl6bTg+kTvj0v7611H6j8UOBe4AvgB8MyRcX8AnN0/dzJwp/mmM8ay+ivg28DNFihzQ98ATwC+0ZfXecBRI+Vu0ZflJb0fTgN+e7H5mTWtGy3/Pmwr4FvAm/v7xwBrR8b/OXB+r/sc4LHA/sA1wLW9L741sixfB3wZ+CVw19Hl29v55b78Lge+Czx2ZFo/pK9/c6yDP+59dWX/e2iv70sj5R/W++Xy/v9hs9az1/bpXwH8J7DzxojZZdwW3Ki/Rta5E4AP9vk8tPfVqfwmft4ObNPLb937dY/+/n19/Kf6578C7LnUsn38AcD3+vL4+973h84zL9+lx/s84/cBPg38rJd9Sh/+x309vKavFx8do9/uAfwaeMACZd5Hjz/gtsBJ/GZ78nFgl5Gyz+vL4gpaHB7ch+9N225cTttmfmCead0VqDmGHwlcAKS/Xws8pr9+CPB12rbip8Df9OEXzIqTB9G2/af0ZfUzWlw9H/j8rOX6J7Ttx8XAG+jbrb5OHTdXe2kHhb8Gru7Te9sc68lOvT/X9X56+cg8PR/4AvBW2vp5LvA/Nzg2NlIA7g9cB2w9x7j3AB/sr4/qK+mTaUft23Hjjds+vfMeQdsjfXMvv1Ay/udez32BXwH36OMf2FeOrXvZs4GXjLRrScm4Dz8FeOP61E8LnqcAtwR2AD4EfKyPuxVtBb5bf38H4J799UHAGlqwbg28Avivheajr0CPmGfeTgVes8jyHE3GjwHu3ZfXfWhB9uQ+7g9pG4Fb0hLoA2k7ZvPOzxzTumGZzhp+NPDVkTas7a/vRtspuOPIenCX+erqy/LHwD17/23DTZPxdcBL+7in0zZUt+njf8j8yXiP3ldbj4w/lJ6MaUf6lwLP7tM+pL+/7Ujbvk/bQG7X379hY8TsMm4LbtRffdhf0RLTk/hN3D8IeHDvlzvTEuThvfxcCfZiYGVfRieMLIOllL0dLTkd1Me9jLZ9OXSeeTmOtuN6KLDXrHHb03YIn9Pb8EDaTundRtpx1KzPvBN4+zzTOhz4/iJ9O5qMV9DOIm1Hi7l/Az7cx926r8N7jcTfPv31h2g7szej7Uw/fJ5pzZeM9+79PVP3aDI+DTikv94BePB8ddES3nXAH9G2HdsxdzL+NO1A5k607eChI+vUcfO1F/jS6HKdYz35QO+zHWjr3xrguSNtu5Z2ELQVbYfgvA2NjY11Sm5n4OKqum6OcT/p42d8pao+VlXXV9UvZ5V9KvDxqvpSVV0DvIrWgQt5TVX9sqq+RTuaui9AVZ1eVadW1XVV9UNaIDx66bN2IxfQNrBLrr+qLqmqj1TVL6rqCtrR2mj564F7Jdmuqn5SVWf24S8E/rqqzu79+3rgfknutMC0dqqq+W5Kui1tmYylqj5fVd/uy+sM2tHNTLuv7fXdtdo1mdOr6ueLzM+4bujrWX4N3BzYJ8k2VfXDqlrs1OFxVXVmX1bXzjH+IuBtVXVtVZ1AO9p+whLbO5cnAP9dVe/t0/4g7QjqSSNl/l9Vfa/HwonA/SYw3Wn0par6+EzcV9VpVfXV3i/nAseycHx+uKpW9+X3fhbup/nKPhH4ZlX9ex/3Vlrins8f05L5i4Gzk/x3kpn7AQ4CvldV/9rn4XTgY7Rt2Jyq6g+r6sXzjF5qXK6rqo/2vvw5bbsw2n9Fi79b9Pg7qw+/lrYTeYequrqqvjzuNLsL+v+5YvNaYK8kt62qK6rqq4vU9eOq+se+7ZidC2a8oaouraof0Y6iD1lie28iyTbA04AjezvPpa0Lzx4p9v2qendV/Zp2QLlrkp3nqG5sGysZXwzsPM95/ztw4xX+vAXquePo+Kr6BW1vcyEXjrz+BW2PlSR7J/nEzM0XtJV1gzoT2IV2SmXJ9Se5ZZJ3JvlRL38KsFOSrarqKtoR2QuBnyT5ZJK794/eCfi7JJcluaxPP70t6+MS2jIZS5IHJ/lcknVJLu9tnJnP99JOmx+f5IIkb+oJcqH5GdcNfT2qqtYAL6EdoV6U5Pgkd1ykroXWOYDzq+8Sdz+irYsb6o69rlE/4sbLbs71dzN0o2WQ5O59vZiJn6NZOD6X0k/zlZ29fSnakd2c+o7zX1XVA2jJ8t+AjyTZkRaXD5+Jyx6bT2cJsTXLUuNy+yTvSvLj3n+fpfdfT86H0K49X9i3U3v3jx5BOyuwOu1m2+cusZ0z6+5NYhP437Szm+ck+VqSxa51LxaXs8tMKi5vRzviHY3NxeISNjA2N1Yy/grtFPH/Gh2YZHvaNZrPjAxe6Ej3J8CuI5/fjhYE6+MfaUche1XVrYG/oCWx9ZJkN9qpqJk7B5da/xG0U6wP7uUfNVM1QFWdXFWPowXkd2mn36GtjH/Yj3Zn/rarqv9az1n5NPB7S7iR6QPAKmC3qtoR+KeRNl9bVa+pqn1o10afSDttt9D8LKq37Un8pq9vpKo+UFWPoG0Qi3aNCOZftxY7u7LLrBvzduc3RwBX0U7Dz7j9Euq9oLdx1O6005tbmtl99U7gO7SzKremnQVb7psjZ29fxt6prarLgb+mbZD3oMXlZ2bF5fZVdfjMR5bYts8AeyS5/5jl/wzYE9i399/vzGrvp6pqP1r8raH1N/0o+flVdQdasj42yZ5LaOfv0RLVmtkjquqcqjqYluz+lrbjcgvWPy4Bdht5PW5cLlb3RbQzbKOxuexxuVGScV9RXwP8fZL9k2yT9v2wE2l7nu8ds6oPA09K8rB+B/ZRrH+A7kC7bnllPyr7o/WppB/RPhr4d+BrtJsmxqn/p7RrEaPt+SVwWb+j+dUj0/jtJAelfbXoV7Tr5tf30f8EvDzJPXvZHZP8/gLTWcxbaNeU3jNzqjvJLknekuQ+c5TfAfhZVV2dZF/gGSPt/h9J7p1kq94X1wLXLzI/80r7buE9aKfCb9/bOrvM3ZL8TtrXzK6m9elM3T+lbdCWut7fDnhxX29/n3Z9fmY5fxM4uI9byY1PQ67r056v/08C9k7yjD5vT6cdOXxiie3bHO1Au655VV/mf7gRpvkJ4AFJntTP4v0p7drrnJK8OsnKJNv2pPJi2hHhf9N2UO/Zl+02/W/fJDNfwVtSXFbV2bRT9SckeXSf5na9/rm+TrUD7Yjt0iS3pe3MzLT7Dn0eb0m7Vn8VPUaSPC3JzA7IZbSk9evF2tdj+sW0e1b+fNaZpJkyz06yc1VdT1u21ad7EVBJlrKdmvF/k+yU9v3gF9MuG0CLy0cn2S3t67NHzvrcvP3fL1F8GHh9P8OwJ+2ekUW/FrkhNtrXOKrqTbSjwzfTNsxfpe09PraqfjVmHWfSLpYfT9uLvZK2IMf6/Cz/h5Y4rqAdlZ2wcPGbeEeSK2gL9W3AR4D9+4o2Tv1H0RLeZUme1uvYjnbK/lTgP0bK3ox2M8kFtGB/ND25V9VHaUd+x/fTUd+hnW2Ybzoz3/uc84v4VfUz2lHstcBX+zx+hhY8N9nbpV03O7qXexVtB2vG7Wkr9c9pN7B9gbbjNe/8zOPpSa7sbVhFO2X3wKq6YI6yN6fdVXkxbQ/9drQ7IaHdnAJwSZKvLzC92b4K7NXrfB3w1KqauTzyStpXvi6l7XB+YOZD/TLK64Av9/5/yGilvY4n0s6KXAL8X9rduQtdp9xSHEH7JsQVtKO2pcbnklXVT2mnkt9CWx53oX1TYKHty3t62QtoNxI+oZ++vhz4XeBZtG3VhbQj55lnEbwLuG+SS5N8GKCfVn7HAtN6Ee2M2z/S1rf/Bg4EPjlH2bfQvkJ6CfBftLvHZ2xFO3L+SR//sF43tJvmTktyFe20+4tqge/e9m3JlcAZfX7/V1X96zzFH0+7tn4FLQ88vaquqXaPzF/TtjeX9Z3acX2clni/AXyUdlMdtO3nR2k32H2Ntt0Y9TbgkD69m+zU07Zr19BuOPwCbTnPN18TkTl2YDYZ/TT3ZbRTwT8Yuj2SNh/9jM4FtJ2vBR9cIW2oTe4BBzOnV/opzjfT9nx+OGyrJG0O+mW0nfpljlfSzhB9beBmaQuwaDJO8u4kFyX5zjzjk+TtSdYkOSPJAybfzBs5iLa3egHt1OHBc12fkHRTUxjP0+YRtIc4rKOddv29cS+jSRti0dPUSR5Fuzb7r1V1k2egpt2e/ie06wEPBv6uqh68DG2VtIGMZ2k6LXpkXFWnMPd3xmYcRAvsqqpTad+NnegzgyVNhvEsTadJXDPehRt/8Xot6//ACUnDMp6lAWzUX8JIchhwGMCtbnWrB9797kt96JK05Tn99NMvrqp5v+86FONZWpqFYnkSyfh8bvwUlF2Z50klVXUs7YvrrFy5slavXj2ByUubtySzH5m5nIxnaZksFMuTOE29CnhOvwvzIcDlVTX2A80lTRXjWRrAokfGST5Ie7LMzknW0h7TuA1AVf0T7ZF+j6c9nekXtIeBS5pCxrM0nRZNxlW14E9S9e/4vmihMpKmg/EsTadN7glckiRtbkzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwMZKxkn2T3JOkjVJjpxj/O5JPpfkG0nOSPL4yTdV0oYylqXptGgyTrIVcAxwALAPcEiSfWYVewVwYlXdHzgY+IdJN1TShjGWpek1zpHxvsCaqjq3qq4BjgcOmlWmgFv31zsCF0yuiZImxFiWptQ4yXgX4LyR92v7sFFHAc9KshY4CfiTuSpKcliS1UlWr1u3bj2aK2kDTCyWwXiWJmlSN3AdAhxXVbsCjwfem+QmdVfVsVW1sqpWrlixYkKTljRBY8UyGM/SJI2TjM8Hdht5v2sfNup5wIkAVfUV4BbAzpNooKSJMZalKTVOMj4N2CvJnkm2pd3UsWpWmR8DjwVIcg9aAHveSpouxrI0pRZNxlV1HXA4cDJwNu1OyzOTHJ3kwF7sCOAFSb4FfBA4tKpquRotaemMZWl6bT1Ooao6iXYzx+iwV428Pgt4+GSbJmnSjGVpOvkELkmSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWBjJeMk+yc5J8maJEfOU+ZpSc5KcmaSD0y2mZImwViWptPWixVIshVwDPA4YC1wWpJVVXXWSJm9gJcDD6+qS5PcbrkaLGn9GMvS9BrnyHhfYE1VnVtV1wDHAwfNKvMC4JiquhSgqi6abDMlTYCxLE2pcZLxLsB5I+/X9mGj9gb2TvLlJKcm2X9SDZQ0McayNKUWPU29hHr2Ah4D7AqckuTeVXXZaKEkhwGHAey+++4TmrSkCRorlsF4liZpnCPj84HdRt7v2oeNWgusqqprq+oHwPdoAX0jVXVsVa2sqpUrVqxY3zZLWj8Ti2UwnqVJGicZnwbslWTPJNsCBwOrZpX5GG1PmiQ70051nTvBdkracMayNKUWTcZVdR1wOHAycDZwYlWdmeToJAf2YicDlyQ5C/gc8GdVdclyNVrS0hnL0vRKVQ0y4ZUrV9bq1asHmba0KUlyelWtHLodCzGepcUtFMs+gUuSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGNlYyT7J/knCRrkhy5QLmnJKkkKyfXREmTYixL02nRZJxkK+AY4ABgH+CQJPvMUW4H4E+Br066kZI2nLEsTa9xjoz3BdZU1blVdQ1wPHDQHOVeC7wRuHqC7ZM0OcayNKXGSca7AOeNvF/bh90gyQOA3arqkxNsm6TJMpalKbXBN3AluRnwFuCIMcoelmR1ktXr1q3b0ElLmqClxHIvbzxLEzJOMj4f2G3k/a592IwdgHsBn0/yQ+AhwKq5bvyoqmOramVVrVyxYsX6t1rS+phYLIPxLE3SOMn4NGCvJHsm2RY4GFg1M7KqLq+qnatqj6raAzgVOLCqVi9LiyWtL2NZmlKLJuOqug44HDgZOBs4sarOTHJ0kgOXu4GSJsNYlqbX1uMUqqqTgJNmDXvVPGUfs+HNkrQcjGVpOvkELkmSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYGM99GMwydAtmLyqJX8kr9n8+qFevR79sPl1w/qsDpuuzW0BrufC29zieX1iGVwdZvPIWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBjZWMk+yf5Jwka5IcOcf4lyU5K8kZST6T5E6Tb6qkDWUsS9Np0WScZCvgGOAAYB/gkCT7zCr2DWBlVd0H+DDwpkk3VNKGMZal6TXOkfG+wJqqOreqrgGOBw4aLVBVn6uqX/S3pwK7TraZkibAWJam1DjJeBfgvJH3a/uw+TwP+NSGNErSsjCWpSm19SQrS/IsYCXw6HnGHwYcBrD77rtPctKSJmixWO5ljGdpQsY5Mj4f2G3k/a592I0k2Q/4S+DAqvrVXBVV1bFVtbKqVq5YsWJ92itp/U0slsF4liZpnGR8GrBXkj2TbAscDKwaLZDk/sA7acF70eSbKWkCjGVpSi2ajKvqOuBw4GTgbODEqjozydFJDuzF/gbYHvhQkm8mWTVPdZIGYixL02usa8ZVdRJw0qxhrxp5vd+E2yVpGRjL0nTyCVySJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwEzGkiQNzGQsSdLATMaSJA1srGScZP8k5yRZk+TIOcbfPMkJffxXk+wx6YZK2nDGsjSdFk3GSbYCjgEOAPYBDkmyz6xizwMuraq7Am8F3jjphkraMMayNL3GOTLeF1hTVedW1TXA8cBBs8ocBLynv/4w8NgkmVwzJU2AsSxNqXGS8S7AeSPv1/Zhc5apquuAy4HbTqKBkibGWJam1NYbc2JJDgMO62+vTHLOxpz+InYGLl72qUz/QcZG6YccNdX9sHHWBcZeHe60zM1YL1Mcz8ZyYyw3G6cfNjCWx0nG5wO7jbzftQ+bq8zaJFsDOwKXzK6oqo4Fjh1jmhtdktVVtXLodgzNftis+2BisQzTG8+b8fJbEvuh2VT6YZzT1KcBeyXZM8m2wMHAqlllVgHP7a+fCny2qmpyzZQ0AcayNKUWPTKuquuSHA6cDGwFvLuqzkxyNLC6qlYB/wK8N8ka4Ge0IJc0RYxlaXrFnd4myWH9tNsWzX6wDzZ1Lr/Gfmg2lX4wGUuSNDAfhylJ0sA2i2Sc5PZJjk/y/SSnJzkpyd5J9kjynWWa5tQ9NnCgfnhUkq8nuS7JU5djGktszxB98LIkZyU5I8lnkkzlV5E2BcbyDW3a4mMZtqx43uSTcX860EeBz1fVXarqgcDLgd9e5klP1WMDB+yHHwOHAh9Y5uksasA++AawsqruQ3tq1ZuWeXqbJWO5MZabLS6eq2qT/gN+BzhlnnF7AN8Zef1F4Ov972F9+B2AU4BvAt8BHkm70/S4/v7bwEvnqPtk4KH99da0L5VnS+uHkWkcBzx1S1wXZk3n/sCXh+yHTfXPWJ6O9XgaYnka+qHXsdHieaM+gWuZ3As4fYxyFwGPq6qrk+wFfBBYCTwDOLmqXtcfpH9L4H7ALlV1L4AkO81R340eG5hk5rGBG+XJTXMYqh+myTT0wfOAT63vDGzhjOVmGtbjaTAN/bDR4nlzSMbj2gZ4R5L7Ab8G9u7DTwPenWQb4GNV9c0k5wJ3TvL3wCeB/xykxcvDflimPkjyLNpG4NHL2nq5Djf2Q7NZxPMmf80YOBN44BjlXgr8FLgvrYO3BaiqU4BH0R4DeFyS51TVpb3c54EXAu+ao74bHi2YRR4buJEM1Q/TZLA+SLIf8JfAgVX1qw2bjS2WsdwYy80WFc+bQzL+LHDztIfWA5DkPkkeOavcjsBPqup64Nm0awf0O+V+WlX/TFswD0iyM3CzqvoI8ArgAXNMd9oeGzhUP0yTQfogyf2Bd9IC96JlmK8thbHcGMvNlhXPG+PC9HL/AXcETgS+T9ub+iSwFze+yL8XcAbwLdrdklf24c+lXcz/Bu0mgD1pe05fp134/yZwwBzTvAXwIWAN8DXgzltoPzyI9lN8V9GOJs7cAvvg07Q985kyq4ZeFzbVP2N50H6YqlgesB8GiWefwCVJ0sA2h9PUkiRt0kzGkiQNzGQsSdLATMaSJA3MZCxJ0sBMxpIkDcxkLEnSwDb5ZNx/17L6Y+xI8qkkzx2n7HpM6y+SbAqPkdsgSY5K8r6h2zFqoeW6HnU9Msk5I+9/2B9/NxFJzkzymEnVpy1LkrsnuW7odoxK8rwkH59gfd9P8tD++g2T3K4meU2Sd0yqvo1l8GSc5D+SHD3H8IOSXLjUxFlVB1TVeybQrsckWTur7tdX1fM3tO45prVtkr9NsjbJlT05vG3Mz65X4kzyjCSr+/R+0pPdI5be+g3Xd5Cu6m25JO0HvZ8+Wmbc5drruutCZarqi1V1tw1td5/ecUn+alb996yqz0+i/i1ZXx9m/q5P8suR98/cgHpP7T8CsFCZP0ryvT6tC5N8Isl2Y9S9f5I1Y5TbJ8m/9fX9siTfTPLiJFnKvExCkuOT/CrJFf3vjCSvTbL9TJmq+peqetKYdb1isXLVfp/4KxNo+036u6peXVWHb2jdG9vgyRh4D/CsOVbCZwPvr6qp2kNcJi+nPeB8X2AH4DG0R7YtiyQvA94GvJ72Q927A/8AHLRc0xzDfatqe+ButN8bfUeSV096Iut7VkQbX1VtP/NH++H7J40Me/9yTTfJ79KeW/yUPu17Af82wfrvDnwFOAe4Z1XtRPu5v0cCN5/UdJbotVW1A7ACeAHwP4AvJrnFJCdi/C1gCp49uh1wOfCokWG/BVxN20ADPIH2fNGf03539KiRsnsABWzd338eeH5/vRXwZtrvkp4LvGhW2f8NnA1c0cf/YR9+K+CXwPXAlf3vjsBRwPtGpn0g7Xmpl/Xp3mNk3A+B/0N7ZurlwAnALebpg08AL1nk+awfAdYBPwBe3IfvD1wDXNvb+K0x+nvHXvb3Fygzez4/BFzY5+MU2gZkZtzjgbN6H54P/J8+fOc+X5cBP6M9G/Zm80yvgLvOGvbUvg7cdo7lelfgC709FwMn9OGn9Lqu6vP4dNqOzVrgz/s8vHdm2Kxl9fI+H5cC/29mWQGHAl+aq73AYb3vr+nT+/hIffv11zen7fhc0P/eBty8j5tp2xG032T9CVpEaVsAABkxSURBVPC/h47Jafwb7dORYVsBr6TF7sXA+4GdRmL4+L7uXQZ8lbZd+Vvaz+xd3ZfZ384xrVcAxy/Qlu36cjyvr1N/35fzbbnpduO2c3z+w8BHFqj/7sB1I+//EPguLcbWAH8wMu72wH/0ebyE9iMXM+Ne2depn9O2c4+cZ3rHA6+YNey3aNubmZh7IfDpkX4/po+/nPZM6LsBL+7x8Ks+7x/q5S+kbQvPBH4xMuwR/fUbaL9B/JE+j6fRtzG054YXsOvs9s7X372+d42Ufwotti+jPXd6r5FxF9J+9ek7fV7eD2w7xDo++JFxVf2S9iDw54wMfhrw3ar6Vn9/VR+/Ey0x/1GSJ49R/QuAJwL3px15PnXW+Iv6+FvTEvNbkzygqq4CDgAuqN/siV8w+sEke9NWoJfQ9iZPAj6eZNtZ87E/7QHl96Ft2OdyKvCyJH+c5N6jZwmS3Az4OG2F3wV4LPCSJL9bVf9BO7o9obfxvv0zRyb5xDzTeihtBf/oPOPn8inaw9hvRztiHz0q+RfaTswOtCOIz/bhR9ASzQra0fdf0IJqXP9O+73tfecY91ra75D+FrArbWNIVT2qj79v748T+vvbA7cB7kRLoHN5JvC7wF1ov4c6zqm2Y2l98aY+vblO4/0l8BDaj5rft8/PaN23p+0g7UL7IfNjkvzWYtMW0Dbw/xN4BG09uBZ4ax/3fNr6swttx/Bw4JqqOoK2sX9+X2ZHzFHvqcCBSV6V5KGzYhrgLX1696Ylob2BI6vqEuD3gHNHthuXJHlskgtHPr8fLSGP6ye07dGtaUnxmCT37OP+nHaEvTNwB9qONEnuS9um3Y+2fj2BFo9jqfZTg5+jHa3P9kTarx3dhRaDzwAuraq30xLqa/u8//7IZ54OPI6WLOfyFNpZ0tvQYv/fkmy1SBvn7O/RMknuTTvT9se07dcXgH+fdYT+VNp29a7Ag/v8bHSDJ+PuPcBTR06JPKcPA6CqPl9V366q66vqDFoSHOcHn58GvK2qzquqnwF/PTqyqj5ZVd+v5gu0DfxcK99cng58sqr+v6q6lnYEvh3wsJEyb6+qC/q0P04LjLn8Ne3XRp4JrAbOH7lZ6UHAiqo6uqquqapzgX8GDp6vYVX1hqp64jyjbwtcXEs4/V9V766qK6r9rudRwH2T7NhHXwvsk+TWVXVpVX19ZPgdgDtV1bXVrtOOnYx7n15MC87ZrqUl1jtW1dVV9aVFqrseeHVV/arv/M3lHSPryeuAQ8Zt6yKeCRxdVRdV1TrgNbRLMDOu7eOvraqTaHv3E7mevQV4IS0JXlBVV9P69ul9Z/Za2o7gXarquqo6re9kL6qqPk2LrwcDJwMXJ3ljkpv1jfjzgD+tqsuq6nLakdhC8fiZqro9QE8wO9IS7FiqalVV/aBvpz5NSygz93dcSztztnvfPpzSh19H2x7tA2xVVedW1Q/GnWZ3AfPH361pR/BU1Zm1+E8NvrUvp/ni77/6fF5L68+dmczPPB4MfLTnkGtoBy8raAdno237aY/Pk5h/O72spiIZ943pxcCTk9yFdvTwgZnxSR6c5HNJ1iW5nBaEO49R9R1pp5Jm/Gh0ZJID+s0cP0tyGe2U6zj1ztR9Q33VfkvzPNqe+IzRveFfANszh6r6dVUdU1UPpx39vw54d5J70JNOv8njst7Ov6Adba6PS4Cdx712k2SrtLsdv5/k57TThfCbfnoKrd9+lOQL6XdIAn9DO6X2n0nOTXLkUhqZZBta0PxsjtH/FwjwtbQ7l/9gkerW9Y31QmavJ3ccu7ELu9F6Mkfdl8zaMZp3PdFv9IS7G3DSSFx8g7ZNuy3tjM0XgA+n3Rj5+sWOtEb1xPAEWjz+PvBHtJ2oOwLbAGeOTPdjtKOucer9Ne106B3GbUuSA5N8bWQ79Tv8Jv5eR0uan0uypt8PQlWdCRzZx1+U5P1JlrrN2IW54+9TtP59J3Bhkn8YvdlrHueNO77HwwVMJgZnb6d/TbuctuTt9HKbimTc/SvtiPhZwMlV9dORcR+g/QD4blW1I/BPtI3xYn5CC9gZu8+8SHJz2imVNwO/Xe0mipNG6l3sKO4CWqKcqW9m43D+GO2aV1X9sqqOoV273Ie2kv6gqnYa+duhqh4/Zjtn+wrtms44p/mhnbI5iHZqbUfaNXro/dSPOA6ibYw+RrvkQD+SPqKq7ky7tv6yJI9dQjsPou3df232iKq6sKpeUFV3pF1P+4csfAf1OH00ez2ZuSxxFXDLmRFJbr/Eum+0nsyqW+upn2U5H/idWbFxi6q6uJ8FeVVV3R14FC2hzhy9LuUMzfVVdTLtfoR70bYp19GOuGemuWNVzZx+HafuT9N2YheV5Fa0ezZeC9yub6c+y2/i7/Kq+tOqulOv8xVJHt7HvaeqHgbcmXZp6q/mmsY8092Jdk/DF2eP60fob6mq+9Muv90X+NOZ0fNUuVi/3BB/fafpjrQ4mbkn5pYjZUdjcKnb6a1oiXiDttPLYdqS8X6067yzv8KyA/Czqro6yb6Mf07/RODFSXbt1+FGj862pd10sQ64LskBtOtPM34K3HbkdOxcdT+hXw/ahnaN9FfAf43ZthskeUnaV6m2S7J1P0W9A21P/2vAFUn+vI/fKsm9kjxopJ179GvLi+qn1V5Fu+705CS3TLJNP0vwpjk+skOfr0toAfH6kXZvm+SZSXbsp5d+TjslTJInJrlr30m5nHbTzPVj9MVt0r62cgzwxpp1DaiX+f0ku/a3l9ICcqbun9I2Pkv1or6e3IZ2nXfmevO3gHsmuV+/jHLUrM8tNr0P0jaQK5LsTOv7qfoO9ybsn4A3JNkNIMntkjypv94v7etDN6Otl9cx5jqS5Kl9HdspzcOAhwOn9vX83cDfJdm5j98tyeNG6r7dIkeKrwQel+R1M0erSe6W5ITc9O7l7WhH4hcB1yc5kJYkZ9p6YJI7z46zPu+P7gcdv+Q3NzotKMkt+jb232mJ7CbrapKHJFnZz65dRUuYGxp/D+vbjG1oZ74uAb7ezzh+G3hm3/YdSLvvZcZi/X0C8HtJHtXrPrLXvXo92rispiYZV9UPaYnsVrSj4FF/DByd5AraxuzEMav9Z9o1n2/Rbjy64esJVXUF7e6/E2kb9GeMTreqvkvbkJ7bT0fd6JRJVZ1DO4r/e9op9ifRvnpxzZhtG/UL2l2eF/a6XkT7WsW5/bTKE2nXMX7Qx7+LdpQKba8Z4JIkX4cbHk7yqfkmVlV/C7yMdiPROtrR9+G0I9vZ/pV2mud82h2Jp84a/2zgh2mnsF9Iu0YK7YavT9OugX4F+Ieq+twCffCtJFfSTm0/H3hpVb1qnrIPAr7ay6+iXb87t487CnhPX2ZPW2B6s32Ads/AucD36UcRVfU94Og+L/8NzL4+/S+0a+aXJZmr//6KFvhn0DYqX2cJRyha0Jtoy+WzfdvwX/zmOuMutIRyBe1O2ZP4zQ7WW4HnJLl0nh3QS2nbnO/TEvm7gddU1Uf6+JfQEtVqWgL8D9rNP9C2Natol20u6zuX+yW5eKbyqjqbdm/JPsDZ/dTz8bSj0F+NNqSqLqbdqPZxWhJ5cp+XGfeg3Wh1Be3o/c3Vvr+7HW2bcjHtaH572k7AfF7Z+/DiPr9fpt19PdflnZ1oN0VdRouXHwF/18cdCzyoz/vxC0xvto8Af0Dr+6fQtn+/7uMOp92jcylt/kdvTr1Jf49WWu0eo+fRTqmvo92odVBN4VdmU+PfUyNJkpbB1BwZS5K0pVo0GSd5d5KLknxnnvFJ8va0O/nOSDKJ29ElLQPjWZpO4xwZH0d7cMV8DqBdH9yL9kCFf9zwZklaJsdhPEtTZ9Fk3L9EPtd3zWYcBPxrv939VGCnJGN/h07SxmM8S9NpEteMd+HGX+hey42/UC1p02E8SwPYqL+gkeQw+rOBb3WrWz3w7ne/+8acvLRJOv300y+uqhVDt2M241lamoVieRLJ+Hxu/PSiXZnn6SbVHqx/LMDKlStr9eqp+961NHWS/GjxUhNjPEvLZKFYnsRp6lW0L9AnyUOAy6tq7IegS5oqxrM0gEWPjJN8kPb4tZ2TrAVeTXs8G1X1T7SnwTye9uSkX9B+tkvSFDKepem0aDKuqgV/Sq4/sP1FE2uRpGVjPEvTySdwSZI0MJOxJEkDMxlLkjQwk7EkSQMzGUuSNDCTsSRJAzMZS5I0MJOxJEkDMxlLkjQwk7EkSQMzGUuSNDCTsSRJAzMZS5I0MJOxJEkDMxlLkjQwk7EkSQMzGUuSNDCTsSRJAzMZS5I0MJOxJEkDMxlLkjQwk7EkSQMzGUuSNDCTsSRJAzMZS5I0MJOxJEkDMxlLkjQwk7EkSQMzGUuSNDCTsSRJAzMZS5I0sLGScZL9k5yTZE2SI+cYv3uSzyX5RpIzkjx+8k2VtKGMZWk6LZqMk2wFHAMcAOwDHJJkn1nFXgGcWFX3Bw4G/mHSDZW0YYxlaXqNc2S8L7Cmqs6tqmuA44GDZpUp4Nb99Y7ABZNroqQJMZalKbX1GGV2Ac4beb8WePCsMkcB/5nkT4BbAftNpHWSJslYlqbUpG7gOgQ4rqp2BR4PvDfJTepOcliS1UlWr1u3bkKTljRBY8UyGM/SJI2TjM8Hdht5v2sfNup5wIkAVfUV4BbAzrMrqqpjq2plVa1csWLF+rVY0vqaWCz38cazNCHjJOPTgL2S7JlkW9pNHatmlfkx8FiAJPegBbC7ytJ0MZalKbVoMq6q64DDgZOBs2l3Wp6Z5OgkB/ZiRwAvSPIt4IPAoVVVy9VoSUtnLEvTa5wbuKiqk4CTZg171cjrs4CHT7ZpkibNWJamk0/gkiRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBjZWMk6yf5JzkqxJcuQ8ZZ6W5KwkZyb5wGSbKWkSjGVpOm29WIEkWwHHAI8D1gKnJVlVVWeNlNkLeDnw8Kq6NMntlqvBktaPsSxNr3GOjPcF1lTVuVV1DXA8cNCsMi8AjqmqSwGq6qLJNlPSBBjL0pQaJxnvApw38n5tHzZqb2DvJF9OcmqS/SfVQEkTYyxLU2rR09RLqGcv4DHArsApSe5dVZeNFkpyGHAYwO677z6hSUuaoLFiGYxnaZLGOTI+H9ht5P2ufdiotcCqqrq2qn4AfI8W0DdSVcdW1cqqWrlixYr1bbOk9TOxWAbjWZqkcZLxacBeSfZMsi1wMLBqVpmP0fakSbIz7VTXuRNsp6QNZyxLU2rRZFxV1wGHAycDZwMnVtWZSY5OcmAvdjJwSZKzgM8Bf1ZVlyxXoyUtnbEsTa9U1SATXrlyZa1evXqQaUubkiSnV9XKoduxEONZWtxCsewTuCRJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGthYyTjJ/knOSbImyZELlHtKkkqycnJNlDQpxrI0nRZNxkm2Ao4BDgD2AQ5Jss8c5XYA/hT46qQbKWnDGcvS9BrnyHhfYE1VnVtV1wDHAwfNUe61wBuBqyfYPkmTYyxLU2qcZLwLcN7I+7V92A2SPADYrao+OcG2SZosY1maUht8A1eSmwFvAY4Yo+xhSVYnWb1u3boNnbSkCVpKLPfyxrM0IeMk4/OB3Ube79qHzdgBuBfw+SQ/BB4CrJrrxo+qOraqVlbVyhUrVqx/qyWtj4nFMhjP0iSNk4xPA/ZKsmeSbYGDgVUzI6vq8qrauar2qKo9gFOBA6tq9bK0WNL6MpalKbVoMq6q64DDgZOBs4ETq+rMJEcnOXC5GyhpMoxlaXptPU6hqjoJOGnWsFfNU/YxG94sScvBWJamk0/gkiRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGNtZDPwaTDN2Cyata8kfyms2vH+rV69EPm183rM/qsOna3Bbgei68zS2e1yeWwdVhNo+MJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGZjCVJGpjJWJKkgZmMJUkamMlYkqSBmYwlSRqYyViSpIGNlYyT7J/knCRrkhw5x/iXJTkryRlJPpPkTpNvqqQNZSxL02nRZJxkK+AY4ABgH+CQJPvMKvYNYGVV3Qf4MPCmSTdU0oYxlqXpNc6R8b7Amqo6t6quAY4HDhotUFWfq6pf9LenArtOtpmSJsBYlqbUOMl4F+C8kfdr+7D5PA/41FwjkhyWZHWS1evWrRu/lZImYWKxDMazNEkTvYErybOAlcDfzDW+qo6tqpVVtXLFihWTnLSkCVoslsF4liZp6zHKnA/sNvJ+1z7sRpLsB/wl8Oiq+tVkmidpgoxlaUqNc2R8GrBXkj2TbAscDKwaLZDk/sA7gQOr6qLJN1PSBBjL0pRaNBlX1XXA4cDJwNnAiVV1ZpKjkxzYi/0NsD3woSTfTLJqnuokDcRYlqbXOKepqaqTgJNmDXvVyOv9JtwuScvAWJamk0/gkiRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgZmMpYkaWAmY0mSBjZWMk6yf5JzkqxJcuQc42+e5IQ+/qtJ9ph0QyVtOGNZmk6LJuMkWwHHAAcA+wCHJNlnVrHnAZdW1V2BtwJvnHRDJW0YY1maXuMcGe8LrKmqc6vqGuB44KBZZQ4C3tNffxh4bJJMrpmSJsBYlqbUOMl4F+C8kfdr+7A5y1TVdcDlwG0n0UBJE2MsS1Nq6405sSSHAYf1t1cmOWdjTn8ROwMXL/tUpv8gY6P0Q46a6n7YOOsCY68Od1rmZqyXKY5nY7kxlpuN0w8bGMvjJOPzgd1G3u/ah81VZm2SrYEdgUtmV1RVxwLHjjHNjS7J6qpaOXQ7hmY/bNZ9MLFYhumN5814+S2J/dBsKv0wzmnq04C9kuyZZFvgYGDVrDKrgOf2108FPltVNblmSpoAY1maUoseGVfVdUkOB04GtgLeXVVnJjkaWF1Vq4B/Ad6bZA3wM1qQS5oixrI0veJOb5PksH7abYtmP9gHmzqXX2M/NJtKP5iMJUkamI/DlCRpYJtFMk5y+yTHJ/l+ktOTnJRk7yR7JPnOMk1z6h4bOFA/PCrJ15Ncl+SpyzGNJbZniD54WZKzkpyR5DNJpvKrSJsCY/mGNm3xsQxbVjxv8sm4Px3oo8Dnq+ouVfVA4OXAby/zpKfqsYED9sOPgUOBDyzzdBY1YB98A1hZVfehPbXqTcs8vc2SsdwYy80WF89VtUn/Ab8DnDLPuD2A74y8/iLw9f73sD78DsApwDeB7wCPpN1pelx//23gpXPUfTLw0P56a9qXyrOl9cPINI4DnrolrguzpnN/4MtD9sOm+mcsT8d6PA2xPA390OvYaPG8UZ/AtUzuBZw+RrmLgMdV1dVJ9gI+CKwEngGcXFWv6w/SvyVwP2CXqroXQJKd5qjvRo8NTDLz2MCN8uSmOQzVD9NkGvrgecCn1ncGtnDGcjMN6/E0mIZ+2GjxvDkk43FtA7wjyf2AXwN79+GnAe9Osg3wsar6ZpJzgTsn+Xvgk8B/DtLi5WE/LFMfJHkWbSPw6GVtvVyHG/uh2SzieZO/ZgycCTxwjHIvBX4K3JfWwdsCVNUpwKNojwE8LslzqurSXu7zwAuBd81R3w2PFswijw3cSIbqh2kyWB8k2Q/4S+DAqvrVhs3GFstYbozlZouK580hGX8WuHnaQ+sBSHKfJI+cVW5H4CdVdT3wbNq1A/qdcj+tqn+mLZgHJNkZuFlVfQR4BfCAOaY7bY8NHKofpskgfZDk/sA7aYF70TLM15bCWG6M5WbLiueNcWF6uf+AOwInAt+n7U19EtiLG1/k3ws4A/gW7W7JK/vw59Iu5n+DdhPAnrQ9p6/TLvx/EzhgjmneAvgQsAb4GnDnLbQfHkT7Kb6raEcTZ26BffBp2p75TJlVQ68Lm+qfsTxoP0xVLA/YD4PEs0/gkiRpYJvDaWpJkjZpJmNJkgZmMpYkaWAmY0mSBmYyliRpYCZjSZIGZjKWJGlgJmNJkgb2/wPM8KKeawxP1gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.subplot(221)\n",
    "plt.bar(np.unique(y), np.unique(y, return_counts=True)[1]/original_size, \n",
    "       color=(\"red\", \"green\", \"blue\"), \n",
    "       tick_label=(\"Class 0\", \"Class 1\", \"Class 2\"))\n",
    "plt.ylim(0,1)\n",
    "plt.title(\"Original Dataset: Class Distribution\")\n",
    "\n",
    "plt.subplot(222)\n",
    "plt.bar(np.unique(y_train), np.unique(y_train, return_counts=True)[1]/train_size, \n",
    "       color=(\"red\", \"green\", \"blue\"), \n",
    "       tick_label=(\"Class 0\", \"Class 1\", \"Class 2\"))\n",
    "plt.ylim(0,1)\n",
    "plt.title(\"Training Set: Class Distribution\")\n",
    "\n",
    "plt.subplot(223)\n",
    "plt.bar(np.unique(y_valid), np.unique(y_valid, return_counts=True)[1]/validation_size, \n",
    "       color=(\"red\", \"green\", \"blue\"), \n",
    "       tick_label=(\"Class 0\", \"Class 1\", \"Class 2\"))\n",
    "plt.ylim(0,1)\n",
    "plt.title(\"Validation Set: Class Distribution\")\n",
    "\n",
    "plt.subplot(224)\n",
    "plt.bar(np.unique(y_valid), np.unique(y_valid, return_counts=True)[1]/test_size, \n",
    "       color=(\"red\", \"green\", \"blue\"), \n",
    "       tick_label=(\"Class 0\", \"Class 1\", \"Class 2\"))\n",
    "plt.ylim(0,1)\n",
    "plt.title(\"Test Set:Class Distribution\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now assess the validation and test accuracy of knn (with `k` varying between 1 and 20) on this dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9b3/8deHQIgQDJIAApFFi1f2LUUtVUCqReoPCyoF6wJW7bV1uXWr1RaV3t7eVtuqvdZbW5XqVSjagt6KWqtwxVpbEAyyuKDGEtawJ+wkn98f50yYDMlkMnNOzjmTz/Px4JHJmZnvfGeGnPf5Lud7RFUxxhhjGtIq6AoYY4wJNwsKY4wxSVlQGGOMScqCwhhjTFIWFMYYY5KyoDDGGJOUBYXxhYiUiciXgq5HJkRkjIiUN+PrrRaRMc31esakyoLCGJeIqIh8zqeyp4vIm8keo6oDVHVxE8sdIyI1IlLl/isXkXki8vkmlHGPiPxPU143Hc31OsZ7FhTGBExEWmdYxEZVzQc6AGcA7wNLRGRcxpUzBgsK0wxEpJ+IfCoi09zfy0TkVhFZKSK7ReT3IpLn3jfGPSq+RUS2isgmEZmRpOzuIvKCiOwQkXUick3cffe4R9dPikil27VT0kA5b7g3S90j86/F3VdvXUSkrYjcLyL/FJEtIvLfInJcip9JmYh8V0RWAntFpHV8d52IjBSRZSKyxy37542VqY5yVZ0J/Bb4SdzrPSgi693y3hGRs9zt44E7ga+577vU3T5DRNa6n9snIvLNuLKKRORPIrLL/dyXiEgr977uIvIHEalwv/Mbk72OiQYLCuMrERkOvALcoKpz4u6aAowH+gCDgelx950IFAA9gG8AD4vICQ28xFygHOgOXAz8h4icE3f/RPcxHYEXgP+qrxBVPdu9OURV81X19ynU5T+BU4GhwOfcx8xsoJ71mQZ8BeioqkcS7nsQeFBVjwdOAeY1oVyAPwLDRaS9+/tSt56dgGeAZ0UkT1VfBv4D+L37voe4j98KXAAcD8wAfuF+lwC34HzmnYGuOAGgblj8L1CK81mMA/5NRL6c5HVMBFhQGD+dhbNzvkJV/5Rw30OqulFVd+DsXIbG3XcYmKWqh1V1IVAF/Eti4SJyEjAK+K6qHlDVd3GOpK+Ie9ibqrpQVauBp4Cm7qDqrYuICHAt8B1V3aGqlTg7wqlNKPshVV2vqvsbeN3PiUiRqlap6ttNrPdGQHACElX9H1XdrqpHVPVnQFvq+UxjVPVFVf3YbaX8H/BnnO8zVrduQC/3c1mizqJxnwc6q+osVT2kqp8Av6Fpn4kJIQsK46d/Bd5qYIB2c9ztfUB+3O/bE46wE++P6Q7EdtIxn+EczTb0OnlNHBNoqC6dgXbAO24XzC7gZXd7qtYnue8bOK2V90VkqYhc0IRywfkMFNgF4Hb1rXW7+nbhtJKKGnqyiJwvIm+7XUu7gAlxj78PWAf82e2WusPd3gvoHvs83OfdidPqMBFmQWH89K9ATxH5hU/lbwQ6iUiHuG09gQ0+vV68bcB+YICqdnT/FbiDyqlqcOlmVf1IVacBXXDGGp6L60ZKxSRguarudccjbsfp7jtBVTsCu3FaHMfUQ0TaAn8A7ge6uo9fGHu8qlaq6i2qejJO197N7sD5euDTuM+jo6p2UNUJjb1fE24WFMZPlTjjEGeLyH96XbiqrgfeAn4sInkiMhjnSDzdKZhbgJNTfO0anG6VX4hIFwAR6SEiX07ztesQkctEpLP7OrvczTWNPEfcOtwNXI1zNA/ObKgjQAXQWkRm4ow9xGwBescGpIFcnK6pCuCIiJwPnBf3OheIyOfc7rfdQLVbt38Ale4g/XEikiMiA+XoVN3E1zERYV+Y8ZWq7gLOBc4XkR/68BLTgN44rYv5wN2q+pc0y7oH+J3bbTIlhcd/F6cL5m0R2QP8hST9/k00HlgtIlU4A9tTGxjLAKe7pwpn/GQpMAgYo6p/du9/Badb7EOcrrkD1O32etb9uV1ElrtdeTfiDKDvBC7FGWuK6YvzXquAvwG/UtVF7jjQBTjjTZ/itLp+i9PNdczrNOXDMMESu3CRMcaYZKxFYYwxJikLCmOMMUlZUBhjjEnKgsIYY0xSmS5GFhpFRUXau3fvoKthjDGR8s4772xT1aQnimZNUPTu3Ztly5YFXQ1jjIkUEfmsscdY15MxxpikLCiMMcYkZUFhjDEmqawZozDGOA4fPkx5eTkHDhwIuiomRPLy8iguLqZNmzZNfq4FhceqDlVx31/v41fLfsX2fdspbFfIt0q+xW2jbiM/tykLixqTnvLycjp06EDv3r1x1u0zLZ2qsn37dsrLy+nTp0+Tn29dTx6qOlTFGb89g5++9VO27duGomzbt42fvvVTzvjtGVQdqgq6iqYFOHDgAIWFhRYSppaIUFhYmHYr04LCQ/f99T4+3vkxB47U/TIOHDnAxzs/5r6/3hdQzUxLYyFhEmXyf8KCwkO/WvarY0Ii5sCRAzyy7JFmrpExxmTOgsJD2/dtT37//uT3G5MNxo4dyyuvvFJn2wMPPMB1112X9Hn5+c4Y3saNG7n44ovrfcyYMWMaPbH2gQceYN++fbW/T5gwgV27diV5RtMMHTqUqVNb1mXALSg8VNiuMPn9xyW/35hsMG3aNObOnVtn29y5c5k2bVpKz+/evTvPPfdc2q+fGBQLFy6kY8eOaZcXb+3atVRXV7NkyRL27t3rSZn1OXLkSOMPakYWFB76Vsm3yGudV+99ea3zuK4k+RGVMdng4osv5sUXX+TQoUMAlJWVsXHjRs466yyqqqoYN24cw4cPZ9CgQTz//PPHPL+srIyBAwcCsH//fqZOnUq/fv2YNGkS+/cfvcjfddddR0lJCQMGDODuu+8G4KGHHmLjxo2MHTuWsWPHAs7yPtu2bQPg5z//OQMHDmTgwIE88MADta/Xr18/rrnmGgYMGMB5551X53XizZkzh8svv5zzzjuvTt3XrVvHl770JYYMGcLw4cP5+OOPAfjJT37CoEGDGDJkCHfccQdQt1W0bds2YmvUzZ49m4kTJ3LOOecwbty4pJ/Vk08+yeDBgxkyZAiXX345lZWV9OnTh8OHDwOwZ8+eOr9nTFWz4t+IESM0aJUHK3XAwwM079/zlHuo/Zf373k64OEBWnmwMugqmhZgzZo1tbfBn3+N+cpXvqILFixQVdUf//jHesstt6iq6uHDh3X37t2qqlpRUaGnnHKK1tTUqKpq+/btVVX1008/1QEDBqiq6s9+9jOdMWOGqqqWlpZqTk6OLl26VFVVt2/frqqqR44c0dGjR2tpaamqqvbq1UsrKipq6xL7fdmyZTpw4ECtqqrSyspK7d+/vy5fvlw//fRTzcnJ0RUrVqiq6iWXXKJPPfVUve/r1FNP1c8++0xfeeUVveCCC2q3jxw5Uv/4xz+qqur+/ft17969unDhQj3zzDN17969deo7evTo2vdQUVGhvXr1UlXVJ554Qnv06FH7uIY+q1WrVmnfvn1r32Ps8dOnT9f58+erquqvf/1rvfnmm4+pf/z/jRhgmTayf7UWhYfyc/N5++q3uWH47bC3M9S0ou2Rztz+hdt5++q37TwK02LEdz/FdzupKnfeeSeDBw/mS1/6Ehs2bGDLli0NlvPGG29w2WWXATB48GAGDx5ce9+8efMYPnw4w4YNY/Xq1axZsyZpnd58800mTZpE+/btyc/PZ/LkySxZsgSAPn36MHToUABGjBhBWVnZMc9ftmwZRUVF9OzZk3HjxrFixQp27NhBZWUlGzZsYNKkSYBzYlu7du34y1/+wowZM2jXrh0AnTp1avRzO/fcc2sf19Bn9frrr3PJJZdQVFRUp9yrr76aJ554AoAnnniCGTNmNPp6qbKg8Fh+bj6X9bgX7tsKs6o5r3Qr946910LCBMKvNkVjLrzwQl577TWWL1/Ovn37GDFiBABPP/00FRUVvPPOO7z77rt07do1rbn9n376Kffffz+vvfYaK1eu5Ctf+UpGZ6K3bdu29nZOTk69YwRz5szh/fffp3fv3pxyyins2bOHP/zhD01+rdatW1NTUwNwTJ3bt29fe7upn9WoUaMoKytj8eLFVFdX13bfecGCwgfl5Udv79kTXD2MCUp+fj5jx47lqquuqjOIvXv3brp06UKbNm1YtGgRn32WfIXrs88+m2eeeQaAVatWsXLlSsDpg2/fvj0FBQVs2bKFl156qfY5HTp0oLKy8piyzjrrLBYsWMC+ffvYu3cv8+fP56yzzkrp/dTU1DBv3jzee+89ysrKKCsr4/nnn2fOnDl06NCB4uJiFixYAMDBgwfZt28f5557Lk888UTtwPqOHTsAZ8zknXfeAUg6aN/QZ3XOOefw7LPPsn379jrlAlxxxRVceumlnrYmwOegEJHxIvKBiKwTkTvqub+XiLwmIitFZLGIFMdtXy4i74rIahH5Vz/r6bX4oNi9O7h6GBOkadOmUVpaWicovv71r7Ns2TIGDRrEk08+yWmnnZa0jOuuu46qqir69evHzJkza1smQ4YMYdiwYZx22mlceumljBo1qvY51157LePHj68dzI4ZPnw406dPZ+TIkZx++ulcffXVDBs2LKX3smTJEnr06EH37t1rt5199tmsWbOGTZs28dRTT/HQQw8xePBgvvCFL7B582bGjx/PxIkTKSkpYejQodx///0A3HrrrTzyyCMMGzasdpC9Pg19VgMGDOCuu+5i9OjRDBkyhJtvvrnOc3bu3JnyDLOUNTaIke4/IAf4GDgZyAVKgf4Jj3kWuNK9fQ7wlHs7F2jr3s4HyoDuyV4vDIPZMT/4wdFG+sknB10b09LUN2BpWoZnn31WL7vssgbvT3cw289FAUcC61T1EwARmQtcCMSPOPUHYnG4CFgAoKqH4h7Tloh1kVnXkzGmud1www289NJLLFy40POy/QyKHsD6uN/LgdMTHlMKTAYeBCYBHUSkUFW3i8hJwIvA54DbVHVj4guIyLXAtQA9e/b0/h2kaX3cu7agMMY0h1/+8pe+lR30kfqtwGgRWQGMBjYA1QCqul5VB+MExZUi0jXxyar6qKqWqGpJ585Jrw3erOJbFIcOwcGDwdXFGGMy5WdQbABOivu92N1WS1U3qupkVR0G3OVu25X4GGAVkNr0hICpHm1R5OY6P61VYYyJMj+DYinQV0T6iEguMBV4If4BIlIkIrE6fA943N1eLCLHubdPAL4IfOBjXT2zZw/s3Qvt20NsgoQFhTEmynwLClU9AlwPvAKsBeap6moRmSUiE92HjQE+EJEPga7Aj9zt/YC/i0gp8H/A/ar6nl919VKs26m4GAoKnNsWFMaYKPN1jEJVF6rqqap6iqr+yN02U1VfcG8/p6p93cdcraoH3e2vqupgVR3i/nzUz3p6KT4ojj/euW1BYcKq6lAVdy+6m873dabVva3ofF9n7l50d8ZXYxSR2qU3wFkNtXPnzlxwwQWZVrleCxYsQER4//33fSm/pQt6MDvrWFCYqPDz0r3t27dn1apVtauwvvrqq/To0cOrqh9jzpw5fPGLX2TOnDm+vQZAdXW1r+WHlQWFx2ID2SedZEFhws3vS/dOmDCBF198EXB25PFnC//jH//gzDPPZNiwYXzhC1/ggw+cIchf/OIXXHXVVQC89957DBw4sM61JepTVVXFm2++yWOPPXbMdTDqW+a7viXBFy9eXKe1c/311zN79mzAWXLju9/9LsOHD+fZZ5/lN7/5DZ///OcZMmQIF110UW39tmzZwqRJkxgyZAhDhgzhrbfeYubMmbXLmQPcddddPPjgg+l8nIGyoPCYtShMVPh96d6pU6cyd+5cDhw4wMqVKzn99KOnUZ122mksWbKEFStWMGvWLO68804AbrrpJtatW8f8+fOZMWMGv/71r2tXX23I888/z/jx4zn11FMpLCysXUfppZde4vnnn+fvf/87paWl3H777YCzzMW3v/1tSktLeeutt+jWrVuj76WwsJDly5czdepUJk+ezNKlSyktLaVfv3489thjANx4442MHj2a0tJSli9fzoABA7jqqqt48sknAWe9qLlz59bpkosKP0+4a5Hig8K9dokFhQklvy/dO3jwYMrKypgzZw4TJkyoc9/u3bu58sor+eijjxCR2gvstGrVitmzZzN48GC++c1v1lnDqSFz5szhpptuApxwmjNnDiNGjKh3me/6lgRPxde+9rXa26tWreL73/8+u3btoqqqii9/+csAvP7667WhkJOTQ0FBAQUFBRQWFrJixQq2bNnCsGHDKCyM3pUuLSg8Vl+LwhYGNGFU2K6QbfsaXpTOi0v3Tpw4kVtvvZXFixfXrnYK8IMf/ICxY8cyf/58ysrKGDNmTO19H330Efn5+WzceMxiDMfYsWMHr7/+Ou+99x4iQnV1NSLCffc1rdssfulvSL789/Tp01mwYAFDhgxh9uzZLF68OGnZV199NbNnz2bz5s213WpRY11PHrOuJxMVzXHp3quuuoq7776bQYMG1dm+e/fu2sHt2FhAbPuNN97IG2+8wfbt2xu9dvZzzz3H5ZdfzmeffUZZWRnr16+nT58+LFmypN5lvhtaErxXr16sWbOGgwcPsmvXLl577bUGX7OyspJu3bpx+PBhnn766drt48aN45FHnO666upqdrtHiJMmTeLll19m6dKlta2PqLGg8FBlpdN6yMuDTp0sKEy43TbqNk454ZRjwiKvdR6nnHAKt426LePXKC4u5sYbbzxm++233873vvc9hg0bVuciQd/5znf49re/zamnnspjjz3GHXfcwdatW1m2bBlXX331MeXMmTOnthsp5qKLLmLOnDkNLvNd35LgJ510ElOmTGHgwIFMmTIl6fLjP/zhDzn99NMZNWpUnWXSH3zwQRYtWsSgQYMYMWJE7RX3cnNzGTt2LFOmTCEnJ6dpH2BIiKZyuaoIKCkp0dgFy4Oydi307w99+8KHH8L8+TB5Mlx4IbgHMMb4bu3atfTr1y+lx1YdquK+v97HI8seYfv+7RQeV8h1Jddx26jb7KqMHqmpqamdMdW3b99A61Lf/w0ReUdVS5I9z8YoPBTf7QTWojDhl5+bz71j7+XesfcGXZWstGbNGi644AImTZoUeEhkwoLCQxYUxph4/fv355NPPgm6GhmzMQoPWVCYsMiWLmXjnUz+T1hQeMiCwoRBXl4e27dvt7AwtVSV7du3p3zeSCLrevKQBYUJg+LiYsrLy6moqAi6KiZE8vLyKI7tnJrIgsJDsXWeYt9Fu3bQqhXs3w+HD0ObNsHVzbQcbdq0oU+fPkFXw2QR63ryUKxFcZJ7XT+Ro62Kyspg6mSMMZmyoPDI3r2wc6dz+dOioqPbrfvJGBN1FhQe2eBeDby42GlJxNh6T8aYqLOg8EjiQHaMtSiMMVFnQeERCwpjTLayoPBIQ0FRUOD8tKAwxkSVBYVH4i+BGs9aFMaYqPM1KERkvIh8ICLrROSOeu7vJSKvichKEVksIsXu9qEi8jcRWe3e97VjSw8X63oyxmQr34JCRHKAh4Hzgf7ANBHpn/Cw+4EnVXUwMAv4sbt9H3CFqg4AxgMPiEhHv+rqBQsKY0y28rNFMRJYp6qfqOohYC5wYcJj+gOvu7cXxe5X1Q9V9SP39kZgK9DZx7pmzILCGJOt/AyKHsD6uN/L3W3xSoHJ7u1JQAcRqXOhXhEZCeQCHye+gIhcKyLLRGRZkOvaHDgA27ZB69bQpUvd+ywojDFRF/Rg9q3AaBFZAYwGNgDVsTtFpBvwFDBDVWsSn6yqj6pqiaqWdO4cXIMjdrJdjx7O2k7xLCiMMVHn56KAG4D4OUDF7rZabrfSZAARyQcuUtVd7u/HAy8Cd6nq2z7WM2OJiwHGs6AwxkSdny2KpUBfEekjIrnAVOCF+AeISJGIxOrwPeBxd3suMB9noPs5H+voicTFAONZUBhjos63oFDVI8D1wCvAWmCeqq4WkVkiMtF92BjgAxH5EOgK/MjdPgU4G5guIu+6/4b6VddMNTSQDRYUxpjo8/V6FKq6EFiYsG1m3O3ngGNaDKr6P8D/+Fk3L6USFLYooDEmqoIezM4K1qIwxmQzCwoPJAuK/Hxn2fG9e6G6+tj7jTEm7CwoPJBs1lOrVtChg3PbrnJnjIkiC4oMHTwIW7dCTg6ceGL9j7HuJ2NMlFlQZGjjRudn9+5OWNTHgsIYE2UWFBlKNj4RY0FhjIkyC4oMWVAYY7KdBUWGLCiMMdnOgiJDFhTGmGxnQZGhZFNjYywojDFRZkGRoWQLAsZYUBhjosyCIkNN6Xqy9Z6MMVFkQZGBw4dh82bn7OuGTrYDa1EYY6LNgiIDmzaBqhMSbdo0/LiCAuenBYUxJoosKDKQSrcTWIvCGBNtFhQZSGXGE1hQGGOizYIiA6nMeAILCmNMtFlQZMC6nowxLYEFRQYsKIwxLYEFRQZSDYr4CxfV1PhbJ2OM8ZqvQSEi40XkAxFZJyJ31HN/LxF5TURWishiESmOu+9lEdklIn/ys46ZSDUocnKgfXtnKu3evf7XyxhjvORbUIhIDvAwcD7QH5gmIv0THnY/8KSqDgZmAT+Ou+8+4HK/6pepI0fqXrSoMdb9ZIyJKj9bFCOBdar6iaoeAuYCFyY8pj/wunt7Ufz9qvoaENqrTG/e7HQjde0KubmNP96CwhgTVX4GRQ9gfdzv5e62eKXAZPf2JKCDiBSm+gIicq2ILBORZRUVFRlVtqlSnRobY0FhjImqoAezbwVGi8gKYDSwAahO9cmq+qiqlqhqSefOnf2qY71SHZ+IsYUBjTFR1drHsjcA8cfbxe62Wqq6EbdFISL5wEWqusvHOnkm3aCwFoUxJmr8bFEsBfqKSB8RyQWmAi/EP0BEikQkVofvAY/7WB9PNTUobGFAY0xU+RYUqnoEuB54BVgLzFPV1SIyS0Qmug8bA3wgIh8CXYEfxZ4vIkuAZ4FxIlIuIl/2q67psBaFMaal8LPrCVVdCCxM2DYz7vZzwHMNPPcsP+uWqVQXBIyxoDDGRFXQg9mRZbOejDEtRaNBISKvpbKtJamubtrJdmBBYYyJrga7nkQkD2gHFInICYC4dx3PsedDtChbtzpnZnfuDHl5qT3HgsIYE1XJxii+Cfwb0B14h6NBsQf4L5/rFWpNHcgGCwpjTHQ1GBSq+iDwoIjcoKq/bMY6hZ4FhTGmJUllMLtGRDrGfhGRE0TkWz7WKfSaOuMJLCiMMdGVSlBcE3+2tKruBK7xr0rhZy0KY0xLkkpQ5IhIbHwitnx4CuulZq+mTo0FW+vJGBNdqZxw9zLwexH5tfv7N91tLVamLQpVOBq9xhgTbqkExXdxwuE69/dXgd/6VqMISCco2rSB446D/fth3z7ninfGGBMFjQaFqtaIyGzgdVX9wP8qhVtNDWxw18Dt0cSzSY4/3gmKPXssKIwx0ZHKmdkTgXdxu5tEZKiIvJD8WdmrogIOHYJOnaBdu6Y91wa0jTFRlMpg9t04lzXdBaCq7wJ9/KxUmKXT7RRjQWGMiaJUguKwqibO1VE/KhMF6cx4irGgMMZEUSqD2atF5FKcabJ9gRuBt/ytVnhZi8IY09Kk0qK4ARgAHASeAXbjrAHVIllQGGNamqQtCvfkulmqeitwV/NUKdwsKIwxLU3SFoWqVgNfbKa6REI66zzFWFAYY6IolTGKFe502GeBvbGNqvpH32oVYtaiMMa0NKkERR6wHTgnbpsCLS4oVC0ojDEtTypjFCtV9RfNVJ9Q274dDh6Ejh0hP7/pz7eFAY0xUZTKGMW0dAsXkfEi8oGIrBORO+q5v5eIvCYiK0VksYgUx913pYh85P67Mt06eCmT1gRAQYHz01oUxpgoSaXr6a8i8l/A76k7RrE82ZPc1sjDwLlAObBURF5Q1TVxD7sfeFJVfyci5wA/Bi4XkU44Z4SX4HRzveM+d2cT3pvnMg0K63oyxkRRKkEx1P05K26bUnfMoj4jgXWq+gmAiMwFLgTig6I/cLN7exGwwL39ZeBVVd3hPvdVYDwwJ4X6Nsnzz8M3vpHaYw8ccH5aUKRm3TqYPBl+8AO45BLvyn32WbjpJmfNLa8NHQovvwytU/nLSMHBg3DeeXD22fDDH3pTJsDq1TB1KvzHf8D/+3/elRu0mhr46lfhrRZ7Sm/6ysshL8+fslNZPXZsmmX3ANbH/V4OnJ7wmFJgMvAgMAnoICKFDTz3mLVaReRa4FqAnj17plXJQ4ecsYdUtWoF5zQWkQ1oaUHx0kvw3nvwzDPeBsXTT8OmTd6VF++11+DDD6F/f2/KW7kS3njDKdPLoPjTn2DVKpg7N7uCorwc/vd/g66FSdRoUIhIAU430Nnupv/DOQnPiyHZW4H/EpHpwBvABqA61Ser6qPAowAlJSVprT81caKzImyq2raFDh3SeaWWFxSxrrrYT6/LffllGDHCu3InT4YlS5zyvQqKWF23bHEOSnI9ujakX59t0GLvZ9gw+POfg61L1LRt61/ZqTSwHwdWAVPc3y8HnsBpCSSzAYhfOq/Y3VZLVTfGyhGRfOAiVd0lIhuAMQnPXZxCXZusbVt/P+B4Le0qd34HxYABUFTkXbm9ex8NCq/EylJ1WkG9enlbbrYGRe/e3n63JjOprPV0iqreraqfuP/uBU5O4XlLgb4i0kdEcoGpQJ3rWIhIkYjE6vA9nFACeAU4T0ROEJETgPPcbZHWtq1zRHn4sNN3ne0Sj6a9cOiQU15ODnTr5k2ZMbEVgf0ICr/KLS93QihbZLI6s/FPKkGxX0Rql/EQkVHA/saepKpHgOtxdvBrgXmqulpEZrkXQwKn1fCBiHwIdAV+5D53B/BDnLBZitPVtSPldxViLan7KfFo2gsbNzo/u3VzwsJLsUkKUQqKQ4dg2zbvyg1apjMLjT9S6Xr6V+BJd6wCYCcwPZXCVXUhsDBh28y4288BzzXw3Mc52sLIGscf7/xh79kDXboEXRv/xJ/FDs5tL7pd/NyRRCUoYq2q+HI7d/am7KBZUIRTKrOeSoEhInK8+3sLOBb2T0tpUWzbVre7yaudpKwxBbsAABjCSURBVAWF06qK724qL3cGf7OBBUU4Ndj1JCI3i0jtGQaqukdV94jIN0SkxV6PIlMtJSjWr0/+e6bl+hkUXtW1pqZuOHhVbmLgeFVuGPj5/Zr0JRuj+DrwZD3bnwKu8qc62a+lrPeUuDOLQouiqMiZbLBrF1RVZV6e360qr8sN2pEjR8eyuncPti6mrmRB0VpVDyduVNVDQJZP7PRPS1nvKbbzir1fr3eSfsyKETkaQBs2JH9sKvz+DLwuN2hbtkB1tTN211xT1k1qkgVFKxHpmrixvm0mdS2l6ym28zr99Lq/e1WuX10TXk6RjZVRUuKE0KZNzlGzV+V6/dkGzabGhleyoLgPeFFERotIB/ffGOBPOIv5mTS0tKA488y6v3tVrl9B4eWAdqyMk0+Grl2dMYvNm70r1+vPNmg2kB1eDc56UtUnRaQCZzHAgTgLAa4GZqrqS81Uv6zT0oJi5Mi6R9OZLLZ3+LBTjoj3J9vF+BEUJ53k/Nu82dmW6Y6wvqDIhjP9LSjCq7HrUbykqqNVtVBVi9zbFhIZaClBEZu94uXR9ObNzg7xxBOhTZvM61gfP4KiuNjbcmOfbb9+zv+n/fthZ6AL8HvDgiK8Ujkz23ioJQRF/Ml2PXp4N+20OaZOejlFNr6+XpWb2KryekpvkGxqbHhZUDSzlhAUO3Y41+4oKHBW2vXqaLo5BjvD3qJIbFX5cZJgUKxFEV4WFM2sJQRF4h+810Hh547Eq1lP8a0qL4PCr882DGzWU3ilHBQicoaIvOxe2/qrflYqm1lQeFeuHzp3do7Ut293+v7TFWtVHX+8P62qbAuKmpqj5670OOYSZSZoyZbwODFh0804V6GbgLOyq0mDBYV35fqhVaujO6pMTrpLPDq2oEhu61ZnVlxRkX+X8zTpS9ai+G8RmSkisa9tF3AxTlhk8W7OXy0hKGKDkrGdpFfdOc3Vh+3FzjexrrHw2bjROfs4XX59tkGz8YlwazAoVPWrwArgTyJyBfBvQFugELCupzS1hKBo6Kg3CrOe4svPpL6JdW3b1lma4siRukuEN1W2tigsKMKtsfMo/hf4MlAAzAc+VNWHVLUJV5k28Y47zjnp7MAB7676FjaJf/SxBd4yOZpuzgXj/GhR+FVufKhF+Up3NjU23JKNUUwUkUXAyzjXzP4acKGIzBWRU5qrgtlGJPtbFYk7s7w8Z5C4ujr9o+nYgnFdu/q/YFyUgqKgANq3h717o70isbUowi1Zi+LfgfOBKcBPVHWXqt4C/AD3kqUmPdkcFInTQmMy3Uk2547Ei35/Pz6D6uqjl4KNtariV7yNcveTTY0Nt2RBsRuYDFwEbI1tVNWPVHWq3xXLZtkcFLt3O0e3+flH3ydEKyi8PPKP3/FlWm5Dy3BnU1BYiyKckgXFJJyB69bApc1TnZYhm4MifgcZv0hdpkfpUQqKhlpVmX4GiTOevCo3DCwowi3Z6rHbgF82Y11ajGwOioYGJaPUoujaFXJynLn9Bw82fUykuVtVUW9RJK4NZsLH1yU8RGS8iHwgIutE5I567u8pIotEZIWIrBSRCe72XBF5QkTeE5FS9zoYWSObg6KxnVm6U06bc1ZMTs7RMYB0TrqLr2t8qyrTz8CvzzZoFRXODMATTnAG5k34+BYUIpIDPIwzIN4fmCYi/RMe9n1gnqoOA6YCv3K3XwOgqoOAc4GfiUjWrEvVkoMi06Pp5hrszKS+DX0G8Wd819R4V27UWxTW7RR+fu58RwLrVPUT9zrbc4ELEx6jQKxxXgC4czroD7wOoKpbcc4KL/Gxrs3KgsK7cv3iR1AcdxwUFjpLhVekcSaSBYUJip9B0QOIbwyXu9vi3QNcJiLlwELgBnd7KTBRRFqLSB9gBHDMsaSIXCsiy0RkWUU6f3kBaYlBkcnRdBALxmUyQJys9eNHAGVLUNjU2PAKujtnGjBbVYtxFht8yu1iehwnWJYBDwBvAcec06uqj6pqiaqWdO7cuRmrnZmWGBTt2kGnTukdTQexYJwfO/RMy21onKZTJ+dz2bMnmv+nrEURfn4GxQbqtgKK3W3xvgHMA1DVvwF5QJGqHlHV76jqUFW9EOgIfOhjXZtVNgdFQ1M447c1dScZxI7Er6BI9zOIb1UllitytNxMVrwNigVF+PkZFEuBviLSR0RycQarX0h4zD+BcQAi0g8nKCpEpJ2ItHe3nwscUdU1Pta1WcWCIspLLtRnzx6orHRaDx07Hnt/urNzglgHKJOZRMnqm265jbWqotz9ZEERfg2eR5EpVT0iItcDrwA5wOOqulpEZgHLVPUF4BbgNyLyHZyB7emqqiLSBXhFRGpwWiGX+1XPIBQUOD+zrUUR/wcfPy00Jt2dWTa1KPz6DKI8RdYWBAw/34ICQFUX4gxSx2+bGXd7DTCqnueVAf/iZ92ClK1dT6nuzNLdSTbnYOeJJzoXMdqyxZnjn5ub2vNirarjjnPOC0jkd1BErUXR0FnsJlyCHsxukSwovC3XD23aOGGhenR581Q0tIRJjAVFXYmXjDXhZEERAAsKb8v1SzoDz035DJpy/YjGumeiGhQ2NTYaLCgCkK1BkWzGU/z2qARFOjvfxuravr3TJXXwIGzf3vRyvf5sg2bdTtFgQRGA9u2dbol9+5yZLNmisT/62MlyTTmaDnLBOD+Cwq9yo96isKAINwuKAMRf5a6yMti6eKmxP/r8fGfa7MGDsG1bamUGuWBcOjOJUpnBk065jX22RUXOgPvOnc7KtVFhQRENFhQBycbuJz+OpoPsww5LiyKVVlVUr3RnU2OjwYIiINkWFFVVsGuXc+2GwsKGH5duUASxIwlLUGzbllqrKopBYS2KaLCgCEi2BUX88hL1TQuNaSlBkawF1NRyUz3qtqAwfrGgCEi2BUWqf/BRCoru3Z3Q27QptUkHQbeqohYU8V1qNj023CwoApJtQdHY1NiYpk7jDDIocnOdy6LW1MDmzY0/PtVWVbqfgdefbdAaumSsCR8LioDE1nvKloUBs7FFEf+6qdQ3nc8glWnC2dqiaGxtMBMeFhQBybYWRVN3ZqlODU21peKXptQ31bGEDh2c73/fPmc6a2NaQlCYcLOgCEhLD4pUjqaDPNkuxo8WhV/lRm0FWZsaGx0WFAFpqUERW/xt//7Gj6ZjC8YVFAS3YFw6O/RUWj9+BEWXLtC6tbM0yP79jZcbNGtRRIcFRUBaalDEP6axnWQYdiRNGSD24zNQTf3Iu1WrutcmDzub8RQdFhQByaag2LfPOYpt0wZSuXR5lIIi6K6npi7DHaVxijB8vyY1FhQByaagiJ8W2iqF/1GpHqWHYUfiV1A09TNI9ag7SlNkw/D9mtRYUAQkm4KiqX/wqQ66hmGws3t35+eGDVBd3fDj0m1VNfYZpPvZWlAYL1lQBMSCwvujaT/k5Tk7/upq57KoDYm1qnr0SK1V5Vf3W1SCYs8e519Dl4w14WJBERALimh0PcW/frL6NjXU4lsUyaYJ+9VaC5qdbBctvgaFiIwXkQ9EZJ2I3FHP/T1FZJGIrBCRlSIywd3eRkR+JyLvichaEfmen/UMggVFdgZFqnUtKHBWgt27N/n/gaZ2v0WlRRGW79akxregEJEc4GHgfKA/ME1E+ic87PvAPFUdBkwFfuVuvwRoq6qDgBHAN0Wkt191DUJ+vvOzstJZSyjK0t2ZJTuajj/ZLuidSSoDxE2ta6rXj8jWrqcwdCua1PnZohgJrFPVT1T1EDAXuDDhMQrElgMrADbGbW8vIq2B44BDQBYcex+Vk3N0umPUr3LX1D/6jh0bP5oO04JxfuzQoWkBlOpne+KJzv+trVudKwmGVVgOAkxq/AyKHkB8T2m5uy3ePcBlIlIOLARucLc/B+wFNgH/BO5X1R2JLyAi14rIMhFZVlFR4XH1/Zct3U9+HE2HqQ/br6BorNx0WlU5OUdnam3cmPyxQbKgiJagB7OnAbNVtRiYADwlIq1wWiPVQHegD3CLiJyc+GRVfVRVS1S1pHMqcxJDJhuC4sAB57rWrVs7S0ikqrFB16AXA4yXygBxOlN5Gys33VZVFLqfLCiixc+g2ADE/5kXu9vifQOYB6CqfwPygCLgUuBlVT2sqluBvwIlPtY1ENkQFLGj1u7dnaPZVDWlRRE0P2Y9pVJuuq2qKMx8CsM5MiZ1fgbFUqCviPQRkVycweoXEh7zT2AcgIj0wwmKCnf7Oe729sAZwPs+1jUQ2RAU6e7QoxQU8esn1TfxINNWVUOfQbo7U2tRGK/5FhSqegS4HngFWIszu2m1iMwSkYnuw24BrhGRUmAOMF1VFWe2VL6IrMYJnCdUdaVfdQ2KBUU0gqJdO+jUCQ4fdgIhUdhaVWEPivhLxhYVBV0bk4rWfhauqgtxBqnjt82Mu70GGFXP86pwpshmtWwICr+OesMUFOB0Ke3Y4dSra9e694UtLMMeFKleMtaER9CD2S1aNgRFuvPhG5saGragSLbzTbeunTo5S1js3l3/FGm/Ptughe27NY2zoAhQNgVFuke9UZj1BMnrm26rqrFpwtnaorCgiB4LigC15KDo1MlZcC+2OFy8PXucI+x27ZyT88IglR16OqHmR1CceKKzMOHmzc64SthYUESPBUWAWnJQxB9NJ16NLUwn28X4sUNvrNx0Wypt2jhhoRrOk+5samz0WFAEKOpBceiQs/R2q1bOjqmpGtpJhvGIs7mDItaqSncZ7jB3P4Xx+zXJWVAEqKDA+RnVoNi40Tlq7dbNOYegqSwoUvsM0mlVWVAYL1lQBCjWoti9O9h6pCvTP/ioBkX8irdhbVVFISjCMlHBNM6CIkBR73rK9A++oWmcYQyK/HxnYP3gQeeSpzGbNmXWqmrsM/D6sw3a/v1Nu2SsCQcLigBlS1BketSbOOU0bFNjY+qrb6Z1begzyNYWRVMvGWvCwb6qAFlQ1C3Hq3L9Ul99M61rURHk5sLOnc5KsTGZzgwKa1CE9bs1yVlQBCh24aI9e5JfNzmsLCgyr2tD04T9aq0FzabGRpMFRYBat3ZOKlOtezQZFZn+0dd3NB2/YFxhoTf19IofQeFXubGLF23aBEeOpF83r4X1IMAkZ0ERsCh3P2X6R9+qVd0lvON/hulku5goBUVurrN4YU2Nc4Z2WFhQRJMFRcCiGhSHDztHqyJHj17TkTg7J8w7kvpmEnlR38RyvVqGO4wzn2xqbDRZUAQsqkGxebPTZXbiic5Ux3QlHk2HeUfSXC0Kr1pVYRzQDvOBgGmYr9ejMI2LalB49QefOOga5sHO+LqqQnW1N60qvz4DCwrjFQuKgFlQ1C0vzDuS4493ZqpVVh4dgK+pcU6286NVlW1BcfAgbN3qXAUw8eJPJtys6ylgFhR1ywtzUEDd+ob9MwjbFNlYl1pTLxlrgmdBEbDYwoBRW+/Jr+6RlhgUXbo4U6W3bYMDB7K3RRH279Y0zIIiYNaiqFte2Hcm8TOJvKpr4jThbA+KME5UMMlZUAQs6kGR6R99165Hj6Z37nQWjMvNDe+CcX60KKD+AMr0s42Fz8aNzsB70MJ+EGAa5mtQiMh4EflARNaJyB313N9TRBaJyAoRWSkiE9ztXxeRd+P+1YjIUD/rGpSoB4WXR9N//7vzM8wLxsX3+3u5eGF95Wb62eblOYF75IgziBw0C4ro8u3PUURygIeB84H+wDQR6Z/wsO8D81R1GDAV+BWAqj6tqkNVdShwOfCpqr7rV12DFMWgqK4+eonNTKaFxsR2HH/7W93fw8ivFkWsjA8/hB07vFuGO0zdTxYU0eXncdtIYJ2qfqKqh4C5wIUJj1HA3VVSANR3hd9p7nOzUhSDYssWJyy6dHHOHs5UbMfx9tt1fw8jv4PC61aVBYXxgp/nUfQA4ifmlQOnJzzmHuDPInID0B74Uj3lfI1jAwYAEbkWuBagZ8+eGVY3GFEMCq9PikvcSYZ5RxKr2z//6ZwXAN62qrz+DMI0RTbMJ1Oa5ILuCZ4GzFbVYmAC8JSI1NZJRE4H9qnqqvqerKqPqmqJqpZ0DuvoZyOiGBReHxnGyolNEQ7zjqRjR2fF3337/GlVef0ZhKVFEX/J2G7dgq2LaTo/g2IDED/MV+xui/cNYB6Aqv4NyAPil0GbCszxsY6Bi3JQeDXNMbGcMAeFSN36elXXxM/A68826KDI9JKxJlh+BsVSoK+I9BGRXJyd/gsJj/knMA5ARPrhBEWF+3srYApZPD4B0Q4Kr496Y8I+zz6+vl7VtWvXumcrZ1uLwsYnos23oFDVI8D1wCvAWpzZTatFZJaITHQfdgtwjYiU4rQcpqvWXuvtbGC9qn7iVx3DIIpXufM7KMK+M4mvn1d1zcmpO9ZhQWHCxNdGoKouBBYmbJsZd3sNMKqB5y4GzvCzfmHQtq3z7+BBZ/mG444LukaN8/qP/sQTnR1ldbXTLdGlizfl+sWPoIiV5fWAb/wZ3zU1wZ2fYkERbUEPZhuOrvcUle4nr//oc3KODnBGYcE4P4PC63LbtYNOnZzB5IoKb8pMh814ijYLihCIjVNEYWHAmpqjq4DGjla9ENuBRGFH4ndQeL0Mdxi6n6xFEW02/yAEYkHx2WdwwgnB1qUxFRXOZVALC73tJrOgOFqW162q4mJYuRLWrIGgTjcqKztaFxM9FhQhEAuK884Lth5N4fXMpFh5YZ/xBP5Mj40v16/P9oorvC03HVH4fs2xLChC4NJL4f33ncXboqBVK5gxw9syL7oI/vIXmDTJ23L90KkTTJni9P/n5XlX7pgxMHw4XHmld2WCU9eXXnJOEgzSGWdYUESVaFTmZDaipKREly1bFnQ1jDEmUkTkHVUtSfYYG8w2xhiTlAWFMcaYpCwojDHGJGVBYYwxJikLCmOMMUlZUBhjjEnKgsIYY0xSFhTGGGOSypoT7kSkAvgsYXMRsC2A6vjN3lf0ZOt7y9b3Bdn73hLfVy9VTXot6awJivqIyLLGzjiMIntf0ZOt7y1b3xdk73tL531Z15MxxpikLCiMMcYkle1B8WjQFfCJva/oydb3lq3vC7L3vTX5fWX1GIUxxpjMZXuLwhhjTIYsKIwxxiSVlUEhIuNF5AMRWScidwRdHy+JSJmIvCci74pIZK/UJCKPi8hWEVkVt62TiLwqIh+5P0N+BfH6NfDe7hGRDe739q6ITAiyjukQkZNEZJGIrBGR1SJyk7s90t9bkvcV6e9MRPJE5B8iUuq+r3vd7X1E5O/u/vH3IpLbaFnZNkYhIjnAh8C5QDmwFJimqmsCrZhHRKQMKFHVSJ8IJCJnA1XAk6o60N32U2CHqv6nG/AnqOp3g6xnOhp4b/cAVap6f5B1y4SIdAO6qepyEekAvAN8FZhOhL+3JO9rChH+zkREgPaqWiUibYA3gZuAm4E/qupcEflvoFRVH0lWVja2KEYC61T1E1U9BMwFLgy4TiaBqr4B7EjYfCHwO/f273D+WCOngfcWeaq6SVWXu7crgbVADyL+vSV5X5Gmjir31zbuPwXOAZ5zt6f0fWVjUPQA1sf9Xk4WfOlxFPiziLwjItcGXRmPdVXVTe7tzUDXICvjg+tFZKXbNRWp7plEItIbGAb8nSz63hLeF0T8OxORHBF5F9gKvAp8DOxS1SPuQ1LaP2ZjUGS7L6rqcOB84NtuN0fWUadPNJv6RR8BTgGGApuAnwVbnfSJSD7wB+DfVHVP/H1R/t7qeV+R/85UtVpVhwLFOL0tp6VTTjYGxQbgpLjfi91tWUFVN7g/twLzcb78bLHF7S+O9RtvDbg+nlHVLe4fbQ3wGyL6vbl93X8AnlbVP7qbI/+91fe+suU7A1DVXcAi4Eygo4i0du9Kaf+YjUGxFOjrjuznAlOBFwKukydEpL072IaItAfOA1Ylf1akvABc6d6+Eng+wLp4KrYjdU0igt+bOzj6GLBWVX8ed1ekv7eG3lfUvzMR6SwiHd3bx+FM8FmLExgXuw9L6fvKullPAO40tgeAHOBxVf1RwFXyhIicjNOKAGgNPBPV9yYic4AxOEsebwHuBhYA84CeOEvGT1HVyA0KN/DexuB0YShQBnwzrl8/EkTki8AS4D2gxt18J05/fmS/tyTvaxoR/s5EZDDOYHUOTqNgnqrOcvcjc4FOwArgMlU9mLSsbAwKY4wx3snGridjjDEesqAwxhiTlAWFMcaYpCwojDHGJGVBYYwxJikLCmN8JCK941eRNSaKLCiMMcYkZUFhTDMRkZNFZIWIfD7ouhjTFK0bf4gxJlMi8i84Z8NOV9XSoOtjTFNYUBjjv8446+lMzpYLaJmWxbqejPHfbuCfwBeDrogx6bAWhTH+O4Sz+ugrIlKlqs8EXSFjmsKCwphmoKp7ReQC4FU3LLJi6XvTMtjqscYYY5KyMQpjjDFJWVAYY4xJyoLCGGNMUhYUxhhjkrKgMMYYk5QFhTHGmKQsKIwxxiT1/wFpg6g9M63WSwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "k_values = np.arange(1, 30, 1)\n",
    "validation_accuracy = []\n",
    "\n",
    "for k in k_values: \n",
    "    knn = KNeighborsClassifier(n_neighbors=k)\n",
    "    knn.fit(X_train, y_train)\n",
    "    validation_accuracy.append(knn.score(X_valid, y_valid))\n",
    "    \n",
    "\n",
    "plt.plot(figsize=(10,5))\n",
    "plt.plot(k_values, validation_accuracy, \"b-\", linewidth=2, label=\"Validation Accuracy\")\n",
    "plt.plot(k_values[np.argmax(validation_accuracy)], np.max(validation_accuracy), \"og\", markersize=8, label=\"Max. Accuracy\") \n",
    "plt.xlabel(\"k\")\n",
    "plt.ylabel(\"% Correct\")\n",
    "plt.title(\"knn on the Iris Dataset\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The validation accuracy is 0.933 on the validation set, with $k=2$. We now retrain the model with $k=2$ on the training+valdation set, and evaluate it on the test set: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8666666666666667"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "knn_selected = KNeighborsClassifier(n_neighbors=2)\n",
    "knn_selected.fit(X_trainvalid, y_trainvalid)\n",
    "knn.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final accuracy, 0.867, is (expectedly) lower on the test set.\n",
    "\n",
    "**Note**: never use the test set for model selection/tuning. This subset has the only purpose to assess model performance, once the best model has been picked. A large discrepancy between validation and test error should lead us to reconsider the data preprocessing/model building stage, and start it from scratch if necessary. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Cross-Validation\n",
    "An alternative to the holdout method (particularly useful when the dataset is of limited size) is provided by **k-fold-cross-validation**. With this method, instead of partitioning the dataset into a training and validation set, we (randomly) split it into $k$ subsets (called folds). (Do not confuse the $k$ of k-nearest-neighbors with the one of cross-validation!) Iteratively, we retrain the model $k$ times, each time using a different fold as validation set and the remaining $k-1$ folds as training set. This will lead to $k$ different perfomance values (accuracy, error rates, RMSEs, $f_1$ scores, etc.); the final *cross-validation* estimate (which is an estimate of generalization performance) is given by a simple average of these $k$ values. There are no strict rules on how to choose $k$; typical values are $k=5$ and $k=10$ (in `sklearn`, the default value is $k=3$).\n",
    "\n",
    "Therefore, cross-validation works as follows: \n",
    "\n",
    "* First, use fold 1 as test set, and from fold 2 to (k-1) as training set\n",
    "* Next, use fold 2 as test set, and folds 1,3,...,(k-1) as training set\n",
    "* Continue until all k folds are used as test set\n",
    "* Compute the $k$ evaluation scores, and average them \n",
    "\n",
    "Exactly as we have seen for the holdout method, cross-validation can be performed with and without stratification (standard CV is applied in regression, while stratified CV in classification). In the following image, an example of stratified cross-validation is represented. As always, it is preferable to keep a separate test set aside, which will play the same role that we have seen in the holdout method.  \n",
    "\n",
    " <img src=\"./img/model_selection/stratified_cv.png\" width=\"700\" height=\"100\"/>\n",
    " \n",
    "Advantages and disadvantages of CV (compared to holdout): \n",
    "\n",
    " * it offers a more realistic estimate of the generalization performance, since the model is tested on $k$ different subsets of the data, instead of just one \n",
    " * all data points appear in the validation set once\n",
    " * it offers an estimate of the variability of the performance measure (it can be the standard deviation of the measures across the $k$ folds) \n",
    " * more data points are used for training\n",
    " * main disadvantage: computational (each model must be retrained $k$ times) \n",
    " \n",
    "`sklearn` provides us with a function that automatically performs cross-validation for us: `cross_val_score()`. This function takes as main arguments:\n",
    " * estimator: the algorithm to evaluate\n",
    " * X, y: the features and output variable\n",
    " * scoring: the performance measure with which we want to evaluate algorithms\n",
    " * cv: if *integer*, the number of folds, which is 5 by default (if the estimator is a classifier, the function automatically performs stratified cross-validation); otherwise, it can be passed in a form of splitter, such as `KFold` or `StratifiedKFold`\n",
    " \n",
    " The function returns a $k$-dimensional vector of scores. The splitters `KFold` (recommended in regression) and `StratifiedKFold` (recommended in classification) take the following arguments: \n",
    "  * n_splits: the number of folds\n",
    "\n",
    "They create the training/validation sets for each fold. \n",
    "\n",
    "We are now going to see how to perform cross-validation, still using the *iris* dataset. We will perform 5-fold-CV, and will evaluate k-nearest-neighbors with $k=10$. The evaluation will occur with the macro-averaged $f_1$ scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross-Validated Macro f1-scores: 0.98 (+/- 0.02)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "from sklearn import metrics\n",
    "\n",
    "# We set apart 10% of the data for the test set\n",
    "X, y = load_iris(return_X_y=True)\n",
    "X_trainvalid, X_test, y_trainvalid, y_test = train_test_split(X, y, \n",
    "                                                              test_size=0.1, \n",
    "                                                              stratify=y,\n",
    "                                                              random_state=1)\n",
    "# Create 5 folds for cross-validation\n",
    "cv_splits = StratifiedKFold(n_splits=5)\n",
    "# Make f1_macro score for the input of cross_val_score\n",
    "f1_macro = metrics.make_scorer(metrics.f1_score, average='macro')\n",
    "\n",
    "# Perform cv and compute the f1 scores\n",
    "f1_scores = cross_val_score(KNeighborsClassifier(n_neighbors=10), \n",
    "                           X_trainvalid, y_trainvalid, \n",
    "                           scoring=f1_macro,\n",
    "                           cv=cv_splits)\n",
    "\n",
    "# Print cross-validation scores (with 1-sd bounds):\n",
    "print(\"Cross-Validated Macro f1-scores: %0.2f (+/- %0.2f)\" % (f1_scores.mean(),\n",
    "                                                            f1_scores.std()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.9628483, 0.9628483, 1.       , 0.9628483, 1.       ])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f1_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test-set f1-score: 0.93\n"
     ]
    }
   ],
   "source": [
    "# Retrain model on whole dataset and report test error: \n",
    "knn_test = KNeighborsClassifier(n_neighbors=10)\n",
    "knn_test.fit(X_trainvalid, y_trainvalid)\n",
    "print(\"Test-set f1-score: %0.2f\" % metrics.f1_score(y_test, knn_test.predict(X_test), \n",
    "                                                    average='macro') )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If no scoring function is passed as argument, the estimator's default score is used (for classifiers this is the accuracy): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross-Validated accuracies: 0.98 (+/- 0.02)\n"
     ]
    }
   ],
   "source": [
    "# Create 5 folds for cross-validation\n",
    "cv_splits = StratifiedKFold(n_splits=5)\n",
    "# Perform cv and compute accuracies\n",
    "scores = cross_val_score(KNeighborsClassifier(n_neighbors=10), \n",
    "                           X_trainvalid, y_trainvalid, \n",
    "                           cv=cv_splits)\n",
    "\n",
    "# Print cross-validation scores (with 1-sd bounds):\n",
    "print(\"Cross-Validated accuracies: %0.2f (+/- %0.2f)\" % (scores.mean(),\n",
    "                                                         scores.std()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If no splitter is given, but an integer, `cross_val_score` automatically performs cross-validation using that integer as number of folds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cross-Validated accuracies: 0.98 (+/- 0.02)\n"
     ]
    }
   ],
   "source": [
    "# Perform 5-fold-cv (without using the KFold splitter) and compute accuracies:\n",
    "scores = cross_val_score(KNeighborsClassifier(n_neighbors=10), \n",
    "                         X_trainvalid, y_trainvalid, \n",
    "                         cv=5)\n",
    "\n",
    "# Print cross-validation scores (with 1-sd bounds):\n",
    "print(\"Cross-Validated accuracies: %0.2f (+/- %0.2f)\" % (scores.mean(),\n",
    "                                                         scores.std()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the results are the same. There exists also a function, `cross_val_predict` (documented [here](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_predict.html#sklearn.model_selection.cross_val_predict)), that returns cross-validation predictions. With k-nearest-neighbors and $k=10$, this function can be called as follows: \n",
    "```\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "cv_predictions = cross_val_predict(KNeighborsClassifier(n_neighbors=10), \n",
    "                                   X_trainvalid, y_trainvalid, \n",
    "                                   cv=5, method=\"predict\")\n",
    "```\n",
    "\n",
    "In this function (like in cross_val_score), `cv` can take a splitter such as KFold and StratifiedKFold. You can also specify `method=\"predict_proba\"`, in which case the estimated probabilities of each class are returned. The latter can be useful, for instance, to evaluate confusion matrices, ROC and Precision/Recall curves, etc.) of binary classifiers. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Leave-one-out-cross-validation**. Another typical choice for $k$ in Cross-Validation is given by setting $k=n$, where $n$ is the number of observations in the dataset (excluding test data). In this case, we have as many folds as observations, which means that at each cross-validation iteration we use $n-1$ training data. While the method can be computationally expensive (as we have to train $n$ models), it can lead to more accurate results, especially with limited sample sizes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Leave-one-out accuracies: 0.98 (+/- 0.15)\n"
     ]
    }
   ],
   "source": [
    "## Example of leave-one-out CV: \n",
    "from sklearn.model_selection import LeaveOneOut\n",
    "loo = LeaveOneOut()\n",
    "loo_scores = cross_val_score(KNeighborsClassifier(n_neighbors=10), \n",
    "                         X_trainvalid, y_trainvalid, \n",
    "                         cv=loo)\n",
    "print(\"Leave-one-out accuracies: %0.2f (+/- %0.2f)\" % (loo_scores.mean(),\n",
    "                                                       loo_scores.std()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Other Cross-Validation strategies**. There are other cross-validation strategies that can be implemented: \n",
    "\n",
    "* shuffle the data first, and then perform the split; this method will first randomly permute the rows of the dataset, and subsequently will create the folds. It can be activated by setting `shuffle=True` in `KFold` or `StratifiedKFold` (remember to set the `random_state` for reproducibility)\n",
    "* shuffle-split cross-validation: sample an arbitrary number of points for training and validation set, for a pre-specified number of iterations, and then train the algorithms; this method is performed with the `ShuffleSplit` function of sklearn\n",
    "* group cross-validationl; if the data are clustered between groups (for examples, students within schools), it is possible to perform group-specific cross-validation (so that some groups are in the training set, and others in the validation set); the function for this type of split in sklearn is `GroupKFold` \n",
    "* nested cross-validation, which consits of performing an outer loop (in which the test set is rotated) and an innter loop (in which the validation set is rotated). It can be called via `cross_val_score(GridSearchCV())`, for instance; it is more useful for generalization performance estimation than model tuning\n",
    "\n",
    "You can find more in the [documentation](https://scikit-learn.org/stable/modules/cross_validation.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Other remarks on Holdout method and Cross-validation**. \n",
    "* Remember to use the validation set(s) for model tuning and comparison, and the test set only once you have picked the best model. Among the operation that we should perform with the validation set/cross-validation results, we find: (a) tuning (b) model comparison (c) ROC/Precision-Recall curve evaluation (d) choosing an operating point in the ROC/PR curve, and so on. \n",
    "* Once the test set has been set apart, it could seem natural to perform data processing on the remaining part of the data (training+validation) before performing holdout/cross-validation. However, this sequence might easily yield too optimistic results (for example, too large accuracies or too low generalization error estimates). The cause of this behaviour is that we use values of the validation set to modify the training set (for instance, when we perform median imputation of missing data). Instead, we should each time pretend to be in a situation where we only know the training data, and the validation set is not yet known to us. The correct way to go is: (a) perform the splits; (b) perform data preprocessing on the subset assigned for training (in Cross-Validation, this must occur for each cross-validation iteration); (c) use the parameters used to pre-process the training data in step (b) to transform the validation set; (d) perform training and predictions on these pre-processed data. Operations of pre-processing also include feature selection and adjusting class imbalance (with oversampling, for instance). Fortunately, `sklearn` provides us with pipelines that allow to directly incorporate such steps in the holdout/cross-validation routines. In example 2 of the notebook \"Cross-Validation in Practice\" we will see an example on how to perform feature selection with cross-validation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Grid and Random Search \n",
    "**Grid Search**. We have seen an example of grid search above, when discussing the holdout method with knn. Grid search is a method to perform model tuning; it consists of the following steps: \n",
    "\n",
    "1. specify a grid of values for the sets of hyperparameters; the algorithms of step 2 will use such values to train the models \n",
    "1. for each value specified in the grid, train the algorithm and evaluate it on the validation set (alternatively, evaluate it with cross-validation) \n",
    "1. choose the hyperparameter value which leads to best results (e.g., highest accuracy, lowest RMSE, and so on)\n",
    "\n",
    "If it is not sure what values to specify for the initial grid, it is possible to search first across more distant values, in order to have an idea about the *order of magnitude* of the hyperparameter; and then refine the search by looking for values closer to the best performing ones in the previous search. This is known as *telescopic search*. \n",
    "\n",
    "With cross-validation, it is also possible to use the *One-Standard-Error-Rule*: it consists of picking the most parsimonious (=simplest) model whose performance is not 1 standard error beyond the performance measure of the best one. For reference, see [page 244 of The Elements of Statistical Learning](https://web.stanford.edu/~hastie/ElemStatLearn/printings/ESLII_print12.pdf). \n",
    "\n",
    "Here is an example on how to perform grid search with 5-fold-Cross-Validation, in order to pick the best `k` for k-nearest-neighbors with `scikit-learn`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best cv accuracy, best k:\n",
      "0.96,1\n"
     ]
    }
   ],
   "source": [
    "# Telescopic search: we first search for the order of magnitude, and then refine the search \n",
    "grid_init = [1, 20, 50, 100]\n",
    "best_accuracy = 0 \n",
    "\n",
    "for k in grid_init: \n",
    "    knn_tmp = KNeighborsClassifier(n_neighbors=k)\n",
    "    scores = cross_val_score(knn_tmp, X_trainvalid, y_trainvalid, cv=5)\n",
    "    score = np.mean(scores)\n",
    "    # Select the best performing model:\n",
    "    if score > best_accuracy: \n",
    "        best_accuracy = score\n",
    "        chosen_k = k \n",
    "print(\"Best cv accuracy, best k:\\n{},{}\".format(np.round(best_accuracy,2), chosen_k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best cv accuracy, best k:\n",
      "0.985,6\n"
     ]
    }
   ],
   "source": [
    "# Now, we look for values between 1 and 10 to refine the search: \n",
    "grid_refined = np.arange(1, 11, 1)\n",
    "best_new_accuracy = 0 \n",
    "all_scores = []\n",
    "\n",
    "for k in grid_refined: \n",
    "    knn_tmp = KNeighborsClassifier(n_neighbors=k)\n",
    "    scores = cross_val_score(knn_tmp, X_trainvalid, y_trainvalid, cv=5)\n",
    "    score = np.mean(scores)\n",
    "    all_scores.append(scores)\n",
    "    # Select the best performing model:\n",
    "    if score > best_new_accuracy: \n",
    "        best_new_accuracy = score\n",
    "        new_chosen_k = k \n",
    "print(\"Best cv accuracy, best k:\\n{},{}\".format(np.round(best_new_accuracy,3), new_chosen_k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAAFNCAYAAABWoDecAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de3hV1Z3/8feXEAgEBOQiQoBgSdHITY14BwREtA5WUQteUHHqdNo6M6JObWu1Q2d+Y1u10860WlvF2moU8VLaglQFBautBG8gaAkIEgQJd0kghOT7+2PvhEMI5GSTnZOEz+t5zpOzL2evtQN8WHuvtdcxd0dEROqvVaorICLSXClARUQiUoCKiESkABURiUgBKiISkQJURCQiBag0KDNbY2ZjU12PhmJmJ5jZrpiO/Tsz+34cx5bGoQCVZs/MPjCzXeGrwsz2JCx/50iO7e6r3b1DQ9VVWpbWqa6AyJFy95Or3pvZq8Dv3P3XqauRHC3UApXYmNlJZvaxmU0Ol9eY2e1m9r6Z7TCzp80sI9w2ysyKzOw2M9tkZhvM7MYGqsd/mtljCcsDzMwTll83s/8wszfM7HMze9HMjq3vvuH2G83sEzPbbGbfCc9pVBJ1PMbMFprZT8zMGuK8JX4KUImFmZ0KzANucff8hE1XAeOB/sAQ4IaEbT2BTkBv4Cbg52bWJTze1Wb2foxVvhq4HjgOyASm1XdfMxsM/AyYRHAO3QnO6bDMrBuwAJjv7re6nq9uNhSgEofzgNnAFHf/Y41tP3P3T919K/AHYFjCtnJguruXu/scYBcwEMDdn3T3ITHW+RF3X+nupcAzNeqV7L5XAi+4+xvuXgbclUS5WcBCgtsO349efUkFBajE4WvAG+7+ai3bNia8LwUSO2i2uPu+w2yP0+Hqley+vYB1VRvcvQTYVke5/0DQF/GrpGsqTYYCVOLwNaCvmf0k1RUJlQDtE5brvKyOaANBixIAM8sEutTxmYcILt//ZGbt69hXmhgFqMThc4L7nCPM7N5UVwZ4FxhpZn3MrDNwZ0zlPAN82czONLM2wPQkPuME/+GsBmZXdapJ86AAlVi4+3bgAuAiM/vBkR7PzK4xsw8ifvxF4HlgKfAWwf3ZBufu7wO3EgTpp8CW8FVWx+ecoNNsE/C8mbWNo37S8EwdfiLxMLNjgO1AP3dfV9f+0vyoBSrSgMxsgpm1N7MOwP3A2wrPliu2ADWzR8MB0csOsf1EM3vTzMrM7PYa28ab2UdmVmhmcd2vEonDZQSX70VANjA5pbWRWMV2CW9mIwjG8T3u7oNq2d4D6Ad8Gdjm7veF69OAvxPcPysCFgOT3X15LBUVEYkothaouy8Eth5m+yZ3X0wweDrRcKAwnMRhL/AUcGlc9RQRiaop3gPtTcJgZIJWaO8U1UVE5JCa9WxMZnYzcDNAZmbmaSeeeGKKayQiLc2SJUs2u3v32rY1xQBdD/RJWM4K1x3E3R8GHgbIy8vzgoKC+GsnIkcVM1t7qG1N8RJ+MZBjZv3DpzkmEdPAZxGRIxFbC9TM8oFRQDczKwLuAdIB3P0hM+sJFADHAJVm9m9ArrvvNLNvEkyFlgY86u5Rn0AREYlNbAHq7ocd/+buG0mYeKHGtjnAnDjqJSLSUJriJbyISLOgABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAFRGJSAEqIhKRAlREJCIFqIhIRApQEZGIFKAiIhEpQEVEIlKAiohEpAAVEYlIASoiEpECVEQkIgWoiEhEClARkYgUoCIiESlARUQiUoCKiESkABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAFRGJSAEqIhKRAlREJKLYAtTMHjWzTWa27BDbzcx+ZmaFZva+mZ2asK3CzN4NX7PjqqOIyJGIswX6GDD+MNsvAnLC183Agwnbdrv7sPA1Ib4qiohEF1uAuvtCYOthdrkUeNwDfwU6m9nxcdVHRKShpfIeaG9gXcJyUbgOIMPMCszsr2b25cavmohI3VqnugKH0M/d15vZCcB8M1vq7qtq7mRmNxNc/tO3b9/GrqOIHOVS2QJdD/RJWM4K1+HuVT9XA68Cp9R2AHd/2N3z3D2ve/fu8dZWRKSGVAbobGBK2Bt/JrDD3TeYWRczawtgZt2Ac4DlKayniEitYruEN7N8YBTQzcyKgHuAdAB3fwiYA1wMFAKlwI3hR08CfmlmlQQBf6+7K0BFpMmJLUDdfXId2x34Ri3r3wAGx1UvEZGGoieRREQiUoCKiESkABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAFRGJSAEqIhKRAlREJCIFqIhIRApQEZGIFKAiIhEpQEVEIlKAiohEpAAVEYlIASoiEpECVEQkIgWoiEhEClARkYgUoCIiESlARUQiUoCKiESkABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAlWbr/vuDV0vV0s+vJVCAiohEpAAVEYkotgA1s0fNbJOZLTvEdjOzn5lZoZm9b2anJmy73sxWhq/r46qjiMiRiLMF+hgw/jDbLwJywtfNwIMAZnYscA9wBjAcuMfMusRYT5Gk6J6k1BRbgLr7QmDrYXa5FHjcA38FOpvZ8cCFwEvuvtXdtwEvcfggFhFJiVTeA+0NrEtYLgrXHWq9iEiT0qw7kczsZjMrMLOC4uLiVFdHRI4yqQzQ9UCfhOWscN2h1h/E3R929zx3z+vevXtsFRU5GqTiHm9zv6+cygCdDUwJe+PPBHa4+wZgHjDOzLqEnUfjwnUiIk1K67gObGb5wCigm5kVEfSspwO4+0PAHOBioBAoBW4Mt201sx8Ai8NDTXf3w3VGiYgkpaq1e9ttDXO82ALU3SfXsd2Bbxxi26PAo3HUS0SkoTTrTiQRkVRSgIqIRKQAFRGJKLZ7oCItyrZtDP7rLDI/3wjH9IQrroAuesL4aKcWqMjhuMP3vge9ejFq9q2cPe8euPVW6NUrWO+e6hpKCqkFKnI4d98NDzwAe/bQpmpdSUnw84EHgp8/+EEqaiZNgFqgIoeybRvcdx+Ulta+vbQ02L59e+PWS5oMtUBFDmXWLLxVGpvoQSEDKGQAFaSF7wo5ng1YWho88wx89auprq2kgAJUJEFlJRQVwcqVUPiHXhSWfp+ddASgA7toRSV/4wwAMilhQMkqBrzenpyx0LcvpKWlsvbS2BSgclQrL4c1a8LALIRVq2DPnmBb1329OKnNH8nZu4wBFNKTjQAU051CBrCSHArTc3nv41PhXmjTBk44AQYMgJwc6N8f2rZN3blJ/BSg0jxFHFZUWgqrVweBuXIlrF0L+/YF23r1gjPO2B+AXciGXo8Bew44Rg+K6UExZ/MmpGWw47f/ROHm/SH8pz8FnfOtWgWt0pyc4DVgAGRmxnt+chgx/E4VoNK8uAc94/fdx6jKNNL3lsLr7eFf/gVuvx2mTwez6t23bw9CrbAwCLj16/eHW3Y2jB4dhNsXvlBbuHUJjvnAA7V3JLVvD9Om0alfZ07rB6edFqzevTtoyVaVuWABvPRSsO344/eHaU4OHHvskZ2fJCHG36kCVJqXwwwr8vsfYNOuTAqvuLO6hbl5c7BL27bB5fUllwTBlZ2d5OX19OnBz/vuY29lGq33ltIqsz1UVMC0afu3J2jXDgYNCl6w/zZBVaC+9RYsXBhsO/bYAwO158/vxn6iYVMNKsahaOYtZCBwXl6eFxQUpLoaEqdt24Lr7PAmZSVGEVnBvciwb3xnWhe49jo6dG1bHUoDBkCfPkfYwbNtG3++Obj8O2diT7jySujcOdKhKiuDlnBVyBcWws6dQFkZmb/7JQMqPiSHleSwkj6sI43K/R/OyIANGyKXfTgNPdVbkyizxt+ZClpRTjoZlO3fp47fqZktcfe82rapBSrNx6xZkJbGHtqykBG8zFh20AmArmzhJFaQk17EgJMz6Hn7tQ17pdulC0vPDIYqnXOEI5ZatQoCvU+f4BaCOxQXw8r751GYtoyVFX15j6EAtGEvJ7CaARSSw0r6t9pEWw2bSlpZ/nOstlwK6cVKcljNCVzMHC5m7v6djmAomgJUmo2ST7Ywv2Q08zmfUtpzIh8ykWf5In+nC+Fg9jKDPedBM7pNaAY9ekCPDu9xTvmvAWcHx+zv6WcAf+JLOEarUqfvzC+Q0yVCx9RRYNeu/bdKCgvhk+eHUrn7nzAqyaKI81hEDisP/FBpKWzcGKk8Bag0edu3B50wi964hLLWPRi2bzEXMZds1h68c/v20LNn41eyIfTsGdS/pIRO7OQ03uY03gZgNxms4gsUth3Eyk6n1K9jKhmp6PVvgDK3bDkwMDdsCNa3bh0MI7swbys5y37FCbuX0a7GaIpqR/B3RvdApckqLoZ58+DNN4P7hqfnljD+jsH0Kvv40B9qzvcIa9yvq1V4fuWZnQ/omEocv3pQx1TPw3QyJ/RQ7w17qK2qkyyuXv+IZboHf7RV57xyZfArq/q1DBiw/5z79YP0dOr1O9U9UGkR1q+HF1+ExYuD21Nnnw0XXgjdumXCpmvqHFYUR3g2ii7JDZuic2fS2T++9KKLDu6YWrEC/va34GOZmfuDJSenRodaKiZLSbLMigr45JP9rcvCwv27HXPMgeNre/cO7i0fpB6/0yjqDFAzuwX4nbtvi1SCSJI+/hjmzoX33guGGF1wAYwZU+PvdoRhRc1KxPM7ZMfUyv0ttvfeC/atemIq5/hdDPjRbE7YW74/yBJVTZZy220N959S1QQttbQIy2jDx6V9WPnD5RS2283qDe3YuzfY1qMHDB26PzC7d69HwzjGvzN1XsKb2X8Ck4C3Cb7obZ43wet+XcI3T+7w0UdBcH74YdAgGDMGzj+/js6RBhxWlKxGHeYTw/nt2HFgoK5/eQX+xpu02ldGP9aSw0r6sZZ0yvd/KKMd3HEHTJx4hCcUevZZ+PGPYc9uAMpJZw3ZrCSHT+hLJa2w1q3J+tJQBkzKqw7MTp0aoOyIv9MjuoR397vM7HsE389+I/B/ZjYTeMTdV9XzFESAIDjffx/mzAkGmXfqFPQhnHdecEuqTg04rKhJiuH8OnWCvLzgBVC69wVWL3ySleEY2vmMZl/NSNhjMDcbPmuYOvB2NuyZCuxvg7VmH9ms4ULmMYBCvrBvNe1O+zZMqjWzoovhd5rUPVB3dzPbCGwE9gFdgFlm9pK7/3vDVEWOBpWVwb3NF1+ETz+Fbt3g2mvhzDPDm/7SaNr37cagzI8ZVLIMgHJa8xnHUZk4TXC79jDxbrj6tIYp9MmP4IMHYHdwP7IVlRzHZ6Szb/8+mZnNZiRFMvdA/xWYAmwGfg3c4e7lZtYKWAk0vwDdti0YlL1xY/AH1dImakjF+dVRZnl50Js+b17weGWvXnDTTUFrqNab/xK/K64IngcPpbOPLNYfuI9nwNcuhoa6O/JPF8F/3ETNCVoOUFERXF43A8m0QI8FLnf3AwbduXulmV0ST7VikjB8grS04CZ5+xY0UUMqzq+OMsu+O52Fi4yXXgruwWVnw1VXwZAhze9X3ZiPODaKmHuom0yZMUomQOcCW6sWzOwY4CR3/5u7r4itZnFIGD5RrZEmamiUDohUnN8hyiyhPfN/tIL5f36T0mFnc+KJMHUqDBzY/IKzRUvFqIYWNJIimV74d4BTq3rew0v3Anc/tRHql7Q6e+EbYEDtkWhKg7AbdEhKjTK304mXGctCRlBGW4a2Xs5FBT+g/9BjGqbMBKmY/KIxNfde/6ZYZpTf6ZEOpLfEYUvhpXvzG4AfTkRRZTF5bKQno5lPJuGlRHP+fpsa57eGfvyNMw7sVfW2cPvbcP7ohilzwTvgkyGc2aaU9rzLMCpIYzhvMZ4X6dV2J7w1CIY2w9/p0SQVoxpawEiKZIJwtZn9C/BguPx1YHV8VYrJxo0H3HP5mP68whhe4gJGsJCxvEzn0p2RJxVIuY0b8ZJS/s4XmctFrOAk0imnHbv371Nm8GFGMIaiIXyYAWW5VA1JaUUlZ/MGFzKPbmwJ9im15vs7FalDMgH6NeBnwF0E/1JeAW6Os1KxSJioAeAqnuFcXudFxvMKY1jA+ZzV5m0uzMyme4qrWl/usHRfLnPTv8vq8iyOYScTeZYRLDxw3sPMTLj+J/DVsxum4F99AO9+f/991to058k9ROqQzED6TQRPIjVvNYZsAPRiA1OZwQRmM48LeWPfeby+9CpOfwTGjw+er23KKithyZLgKZ71ay+ha+UqruZJzuaNA8fVVWno4SG1/E5jL1OkCUlmHGgGcBNwMlD9jIi7T42xXg3vMMMnurGFa9q/wJe+nssrZ7bltdeCr10YOjSYqKF//xTV+RD27ds/prK4OJjO7MavteX03rtI+58CKK0lPDUkRaTBO+SSuYT/LfAhcCEwHbgGaF7Dl6okDJ84YMxiOHyi8/TvMNGC1ueCBfDKK8EEDCeeGARpqofglJXBokXBPJDbtwdTdn3tazBsWFivM/4DWvkhzy/uISmNVqZIE5FMgA5w9yvN7FJ3/42ZPQksirtisTALxkFOm3bgUzM1hk9kZgZfPjZ27P7A+slPgkHgF1/c+IPAS0qCQJ8/P3g/cCDccEMQ7AfUI8nza1CpKFOkiUgmQKumZtluZoMInofvEV+VGkGXLkkNVcrICKZUGzUK/vrX4PntX/wiGPo4fjycfnq8jyHu2AEvvwyvvRa0PocMCVrCJ5xQxweTPL8GlYoyRVIsmQB92My6EPTCzwY6AN+LtVZNTHp6MEvQOedAQUHQafPoozB7djDR71lnNexEGJs3w5//DH/5S3AlnJcXBHZWVsOVISJH7rABGj51tDOcTHkhUFfbp+bnxwM/BdKAX7v7vTW29yOYY7Q7weOi17p7UbitAlga7vqJu0+oT9lxaNUKhg8PWp5VU7E98QT84Q9BS3XEiCSnYjuEDRuCVu5bbwVXxmefDePGBZPJikjTc9gADZ86+ndgZn0PbGZpwM+BC4AiYLGZzXb35Qm73Qc8Ht5bHQ38N3BduG23uw+rb7mNwSzooR8yZP9kwM8+G/wcPTp4HTAZcB1fnrVmTRCc77wTzBY+enRw/7UlTRAl0hIlcwn/spndDjwNVI+Ydveth/4IAMOBQndfDWBmTwGXAokBmgtMC98vAF5Ist5NglnQkXPiiUEIzp0Lf/xj0Ok0YgSMHeN0fiCYqWhU+OVZvB7MVOS33c7K66YzZ66xYkXQcf2lLwXh2aFDqs9MRJKRTIB+Jfz5jYR1Tt2X872BdQnLRcAZNfZ5D7ic4DL/MqCjmXV19y1AhpkVEEzgfK+7N+lwzc6Gf/7nYJLgF18MhkAtuPctzl66ngv3ZlY/2uglJSxlMHN/BKvnFHDMmNO5/HIYOfLILv9FpPEl8yRSnMPIbyf4ipAbCO6xrgcqwm393H29mZ0AzDezpTW/QsTMbiZ8rLRv374xVjN5vXoF07ZNGLGdebkP8UZ5HosYzuksZiAfsYDzKSKLruVbuPqD73D2vGdI767hPiLNUTJPIk2pbb27P17HR9cDfRKWs8J1icf4lKAFipl1ACa6+/Zw2/rw52ozexU4BVhV4/MPAw9DMJ1dXefSmLoteIZr2jzDl8qf5xXG8BojeYvhHM8GbmQGp7OYtPR28EIznf1JRJK6hD894X0GMIbgGzrrCtDFQI6Z9ScIzknA1Yk7mFk3YKu7VwLfJuiRJxw2VeruZeE+5wA/SqKuTUc4+1NnnIk8x3he5DOOoz8fUz32vbRUMxUdgZY6D2iVln5+LUEyl/C3JC6bWWfgqSQ+t8/MvgnMIxjG9Ki7f2Bm0wkmZJ4NjAL+28yc4BK+6j7rScAvzawSaEVwD3T5QYU0ZTVmf8qklBP4+MB9NFORSLMWZWLkEiCp+6LuPgeYU2Pd3QnvZwGzavncG8DgCHVrOjRTkUiLl8w90D+w/0ucWxEMPar3uNCjjmYqEmnxkmmB3pfwfh+wtuppIalDC/ryLBE5WDJfKtcf2ODue8LldsBx7r4m/uolr84vlUulVHxhl4g0iCP9UrlngMTvgKgI151e++5ykBbw5VkicrBkJmNr7e57qxbC923iq5KISPOQTIAWm1n1TEhmdimwOb4qiYg0D8l+K+cTZvZ/4XIRUOvTSSIiR5NkBtKvAs4MH7XE3XfFXisRkWagzkt4M/t/ZtbZ3Xe5+y4z62Jm/9kYlRMRacqSuQd6UdUEHwDh7PQXx1clEZHmIZkATTOztlUL4TjQtofZX0TkqJBMJ9ITwCtmNgMw4AbgN3FWSkSkOUimE+mHZvYeMJbgmfh5QL+4KyYi0tQlOxvTZwTheSXwMfBsbDVqoTS3o0jLc8gANbMvApPD12aCL5Uzdz+/keomItKkHa4F+iGwCLjE3QsBzOzWRqmViEgzcLhe+MuBDcACM/uVmY2B/d9GISJytDtkgLr7C+4+CTiR4Dvb/w3oYWYPmtm4xqqgiEhTVec4UHcvcfcn3f0fCL5Z8x3gW7HXTESkiUtmIH01d9/m7g+7+5i4KiQi0lzUK0BFRGQ/BaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAFRGJSAEqIhKRAlREJCIFqIhIRApQEZGIFKAiIhEpQEVEIlKAiohEpAAVEYlIASoiElGsAWpm483sIzMrNLM7a9nez8xeMbP3zexVM8tK2Ha9ma0MX9fHWU8RkShiC1AzSwN+DlwE5AKTzSy3xm73AY+7+xBgOvDf4WePBe4BzgCGA/eYWZe46ioiEkWcLdDhQKG7r3b3vcBTwKU19skF5ofvFyRsvxB4yd23uvs24CVgfIx1FRGptzgDtDewLmG5KFyX6D2Cr08GuAzoaGZdk/ysiEhKpboT6XZgpJm9A4wE1gMVyX7YzG42swIzKyguLo6rjiIitYozQNcDfRKWs8J11dz9U3e/3N1PAb4brtuezGfDfR929zx3z+vevXtD119E5LDiDNDFQI6Z9TezNsAkYHbiDmbWzcyq6vBt4NHw/TxgnJl1CTuPxoXrRESajNgC1N33Ad8kCL4VwEx3/8DMppvZhHC3UcBHZvZ34Djgv8LPbgV+QBDCi4Hp4ToRkSbD3D3VdWgQeXl5XlBQkOpqiEgLY2ZL3D2vtm2p7kQSEWm2FKAiIhEpQEVEIlKAiohEpAAVEYlIASoiEpECVEQkIgWoiEhEClARkYgUoCIiESlARUQiUoCKiESkABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCJSgIqIRKQAFRGJSAEqIhKRAlREJCIFqIhIRApQEZGIFKAiIhEpQEVEIlKAiohEpAAVEYlIASoiEpECVEQkIgWoiEhEClARkYgUoCIiESlARUQiUoCKiESkABURiSjWADWz8Wb2kZkVmtmdtWzva2YLzOwdM3vfzC4O12eb2W4zezd8PRRnPUVEomgd14HNLA34OXABUAQsNrPZ7r48Ybe7gJnu/qCZ5QJzgOxw2yp3HxZX/UREjlScLdDhQKG7r3b3vcBTwKU19nHgmPB9J+DTGOsjItKg4gzQ3sC6hOWicF2i7wPXmlkRQevzloRt/cNL+9fM7LzaCjCzm82swMwKiouLG7DqIiJ1S3Un0mTgMXfPAi4GfmtmrYANQF93PwWYBjxpZsfU/LC7P+zuee6e171790atuIhInAG6HuiTsJwVrkt0EzATwN3fBDKAbu5e5u5bwvVLgFXAF2Osq4hIvcUZoIuBHDPrb2ZtgEnA7Br7fAKMATCzkwgCtNjMuoedUJjZCUAOsDrGuoqI1FtsvfDuvs/MvgnMA9KAR939AzObDhS4+2zgNuBXZnYrQYfSDe7uZjYCmG5m5UAl8DV33xpXXUVEojB3T3UdGkReXp4XFBSkuhoi0sKY2RJ3z6ttW6o7kUREmi0FqIhIRLHdA20KysvLKSoqYs+ePamuitQhIyODrKws0tPTU10VkaS16AAtKiqiY8eOZGdnY2apro4cgruzZcsWioqK6N+/f6qrI5K0Fn0Jv2fPHrp27arwbOLMjK5du+pKQZqdFt0CBeoXntu2waxZsHEj9OwJV1wBXbrEVzmppv/kpDlq0S3QpLnD974HvXrBrbfCPfcEP3v1CtYf4VCvF154ATPjww8/bKAK127YsGFMmjQp1jJEZD8FKMDdd8MDD8CePVBSEgRmSUmw/MADwfYjkJ+fz7nnnkt+fn4DVfhgK1asoKKigkWLFlFSUhJbOfv27Yvt2CLNjQJ02za47z4oLa19e2lpsH379kiH37VrF6+//jqPPPIITz31VPX6SZMm8ac//al6+YYbbmDWrFmUlpZy1VVXkZuby2WXXcYZZ5xBMg8I5Ofnc9111zFu3Dh+//vfV68vLCxk7NixDB06lFNPPZVVq1YB8MMf/pDBgwczdOhQ7rwzmOt61KhR1WVt3ryZ7OxsAB577DEmTJjA6NGjGTNmDLt27WLMmDGceuqpDB48+IDyHn/8cYYMGcLQoUO57rrr+Pzzz+nfvz/l5eUA7Ny584BlkWbN3VvE67TTTvOali9fftC6gzz8sHtmpnvQ7qz9lZkZ7BfB7373O586daq7u5911lleUFDg7u7PPfecT5kyxd3dy8rKPCsry0tLS/3HP/6x33zzze7uvnTpUk9LS/PFixfXWc4Xv/hFX7t2rc+bN88vueSS6vXDhw/35557zt3dd+/e7SUlJT5nzhw/66yzvKSkxN3dt2zZ4u7uI0eOrC6ruLjY+/Xr5+7uM2bM8N69e1fvV15e7jt27Kje7wtf+IJXVlb6smXLPCcnx4uLiw847g033ODPP/+8u7v/8pe/9GnTptV6Dkn9eYk0MoJHz2vNHbVAN248dOuzSmlpsF8E+fn51fclJ02aVH0Zf9FFF7FgwQLKysqYO3cuI0aMoF27drz++uvV+w8aNIghQ4bUWUZBQQHdunWjb9++jBkzhnfeeYetW7fy+eefs379ei677DIgGGvZvn17Xn75ZW688Ubat28PwLHHHltnGRWUZ+MAAA0kSURBVBdccEH1fu7Od77zHYYMGcLYsWNZv349n332GfPnz+fKK6+kW7duBxz3H//xH5kxYwYAM2bM4MYbb0z69yfSlLX4Xvg69ewJ7dsH9zwPpX37YL962rp1K/Pnz2fp0qWYGRUVFZgZP/7xj8nIyGDUqFHMmzePp59++og6f/Lz8/nwww+rL7l37tzJs88+W+9jtm7dmsrKSoCDhhRlZmZWv3/iiScoLi5myZIlpKenk52dfdghSOeccw5r1qzh1VdfpaKigkGDBtWrXiJNlVqgV1wBFRWH36eiAq68st6HnjVrFtdddx1r165lzZo1rFu3jv79+7No0SIAvvKVrzBjxgwWLVrE+PHjgSBsZs6cCcDy5ctZunTpYcuorKxk5syZLF26lDVr1rBmzRp+//vfk5+fT8eOHcnKyuKFF14AoKysjNLSUi644AJmzJhBadjy3ro1mOgqOzubJUuWVNf9UHbs2EGPHj1IT09nwYIFrF27FoDRo0fzzDPPsGXLlgOOCzBlyhSuvvpqtT6lRVGAdukCt98etDJr0759sL1z53ofOj8/v/ryucrEiROrL+PHjRvHa6+9xtixY2nTpg0AX//61ykuLiY3N5e77rqLk08+mU6dOgHBpXDNDqVFixbRu3dvevXqVb1uxIgRLF++nA0bNvDb3/6Wn/3sZwwZMoSzzz6bjRs3Mn78eCZMmEBeXh7Dhg3jvvvuA+D222/nwQcf5JRTTmHz5s2HPK9rrrmGgoICBg8ezOOPP86JJ54IwMknn8x3v/tdRo4cydChQ5k2bdoBn9m2bRuTJ0+u9+9RpMk61M3R5vaK3Ink7l5Z6X7XXe4ZGUGHkVnwMyMjWF9ZmdxxGsC+fft89+7d7u5eWFjo2dnZXlZW1mjlx+WZZ57xa6+99rD7qBNJmiIO04mke6AAZvCDH8C0aQc+iXTllZFankeitLSU888/n/LyctydX/ziF9Wt0+bqlltuYe7cucyZMyfVVRFpUArQRF26wFe/mtIqdOzYMalxn83J//7v/6a6CiKx0D1QEZGIFKAiIhEpQEVEIlKA1nD//cFLRKQuCtCYTZ06lR49eiT99E1paSnXXHMNgwcPZtCgQZx77rns2rWL7du384tf/OKQn6uajORwHnvsMbp3786wYcOqX8uXL6/X+YjIfgrQmN1www28+OKLSe//05/+lOOOO46lS5eybNkyHnnkEdLT0+sM0GR95Stf4d13361+5ebmHrC95nR17l79eGddKup6okukhVGAxmzEiBFJTdZRZcOGDfTu3bt6eeDAgbRt25Y777yTVatWMWzYMO644w7cnW9+85sMHDiQsWPHsmnTpsh1fPXVVznvvPOYMGECubm5rFmzhoEDBzJlyhQGDRrEunXryM/Pr24Vf+tb36r+bIcOHbjtttsYOnQob775ZuQ6iDRHR8040JkzYd26uvd79dXkj9mnD1x1VeQq1Wrq1KmMGzeOWbNmMWbMGK6//npycnK49957WbZsGe+++y4Azz33HB999BHLly/ns88+Izc3l6lTp9Z5/KeffprXX3+9erkq9N5++22WLVtG//79WbNmDStXruQ3v/kNZ555Jp9++inf+ta3WLJkCV26dGHcuHG88MILfPnLX6akpIQzzjiD+3XjWI5CaoE2McOGDWP16tXccccdbN26ldNPP50VK1YctN/ChQuZPHkyaWlp9OrVi9GjRyd1/JqX8O3atQNg+PDhB3wjZr9+/TjzzDMBWLx4MaNGjaJ79+60bt2aa665hoULFwKQlpbGxIkTj/S0RZqlo6YFWt+W4m23xVOPZHTo0IHLL7+cyy+/nFatWjFnzpzYQypxurralg8lIyODtLS0OKok0uSpBdrE/OUvf2Hbtm0A7N27l+XLl9OvXz86duzI559/Xr3fiBEjePrpp6moqGDDhg0sWLAgtjoNHz6c1157jc2bN1NRUUF+fj4jR46MrTyR5kIBGrPJkydz1lln8dFHH5GVlcUjjzwCwEMPPcRDDz100P6rVq1i5MiRDB48mFNOOYW8vDwmTpxI165dOeeccxg0aBB33HEHl112GTk5OeTm5jJlyhTOOuus6mPcfffdzJ49u9b6PP300wcMY3rjjTfqPIfjjz+ee++9l/PPP5+hQ4dy2mmncemll0b8jYi0HOZH+JW9TUVeXp7XnIRjxYoVnHTSSfU6TlVfSCov4Y9WUf68ROJmZkvcPa+2bUfNPdBkKThFJFm6hBcRiUgBKiISUYsP0JZyj7el05+TNEctOkAzMjLYsmWL/nE2ce7Oli1byMjISHVVROqlRXciZWVlUVRURHFxcaqrInXIyMggKysr1dUQqZdYA9TMxgM/BdKAX7v7vTW29wV+A3QO97nT3eeE274N3ARUAP/i7vPqW356evoBjyeKiDSk2ALUzNKAnwMXAEXAYjOb7e6JE1DeBcx09wfNLBeYA2SH7ycBJwO9gJfN7IvurvnSRKTJiPMe6HCg0N1Xu/te4Cmg5uMrDhwTvu8EfBq+vxR4yt3L3P1joDA8nohIkxFngPYGEieQKwrXJfo+cK2ZFRG0Pm+px2dFRFIq1Z1Ik4HH3P1+MzsL+K2ZJffdF4CZ3QzcHC7uMrOP6ll+N2BzPT9zJFp6eakoU+U17/JSUWZ9y+t3qA1xBuh6oE/Ccla4LtFNwHgAd3/TzDIITi6Zz+LuDwMPR62gmRUc6hnXOLT08lJRpspr3uWlosyGLC/OS/jFQI6Z9TezNgSdQjWnCPoEGANgZicBGUBxuN8kM2trZv2BHOCtGOsqIlJvsbVA3X2fmX0TmEcwROlRd//AzKYDBe4+G7gN+JWZ3UrQoXSDB6PePzCzmcByYB/wDfXAi0hTE+s90HBM55wa6+5OeL8cOOcQn/0v4L/irB9HcPmv8ppMmSqveZeXijIbrLwWMx+oiEhja9HPwouIxOmoDFAze9TMNpnZskYqr4+ZLTCz5Wb2gZn9a8zlZZjZW2b2Xljef8RZXkK5aWb2jpn9sRHKWmNmS83sXTMrqPsTDVJmZzObZWYfmtmKcOhdXGUNDM+t6rXTzP4trvLCMm8N/74sM7P8cFRMnOX9a1jWB3GdW23/1s3sWDN7ycxWhj+7RC7A3Y+6FzACOBVY1kjlHQ+cGr7vCPwdyI2xPAM6hO/Tgb8BZzbCeU4DngT+2AhlrQG6NfLfm98A/xi+bwN0bqRy04CNQL8Yy+gNfAy0C5dnEnTqxlXeIGAZ0J6gL+ZlYEAM5Rz0bx34EcG8GwB3Aj+MevyjsgXq7guBrY1Y3gZ3fzt8/zmwghifrPLArnAxPXzFerPbzLKALwG/jrOcVDGzTgT/GB8BcPe97r69kYofA6xy97Uxl9MaaGdmrQmC7dM69j8SJwF/c/dSd98HvAZc3tCFHOLf+qUE/xkS/vxy1OMflQGaSmaWDZxC0CqMs5w0M3sX2AS85O6xlgf8D/DvQGXM5VRx4M9mtiR8Ii1u/QnGKM8Ib1P82swyG6FcCMZQ58dZgLuvB+4jGJu9Adjh7n+OschlwHlm1tXM2gMXc+DDM3E6zt03hO83AsdFPZACtBGZWQfgWeDf3H1nnGW5e4W7DyN4imt4fR6RrS8zuwTY5O5L4iqjFue6+6nARcA3zGxEzOW1JrgUfNDdTwFKCC7/YhU+hDIBeCbmcroQtMz6E8yAlmlm18ZVnruvAH4I/Bl4EXiXYOrKRuXBdXzkqzMFaCMxs3SC8HzC3Z9rrHLDy8wFhI/MxuQcYIKZrSGYdWu0mf0uxvKqWky4+ybgeeKfrasIKEpoyc8iCNS4XQS87e6fxVzOWOBjdy9293LgOeDsOAt090fc/TR3HwFsI+gbaAyfmdnxAOHPTVEPpABtBGZmBPfOVrj7A41QXncz6xy+b0cwJ+uHcZXn7t929yx3zya43Jzv7rG1Xsws08w6Vr0HxhFcEsbG3TcC68xsYLhqDMGTcnGbTMyX76FPgDPNrH3493UMwb362JhZj/BnX4L7n0/GWV6C2cD14fvrgd9HPVCqZ2NKCTPLB0YB3cKp9O5x90diLPIc4DpgaXhfEuA7Hs6+H4Pjgd+Ek1q3Ipi0OvahRY3oOOD54N85rYEn3f3FRij3FuCJ8LJ6NXBjnIWF/zlcAPxTnOUAuPvfzGwW8DbB49PvEP8TQs+aWVegnOBx7QbvlKvt3zpwLzDTzG4C1gJXRT5+2JUvIiL1pEt4EZGIFKAiIhEpQEVEIlKAiohEpAAVEYlIASpHFTPLbqxZuKTlU4CKiESkAJWjlpmdEE4Mcnqq6yLN01H5JJJI+EjmUwRzXr6X6vpI86QAlaNRd4Lnny/34IsNRSLRJbwcjXYQTJ5xbqorIs2bWqByNNoLXAbMM7Nd7t5YswBJC6MAlaOSu5eEE0G/FIbo7FTXSZofzcYkIhKR7oGKiESkABURiUgBKiISkQJURCQiBaiISEQKUBGRiBSgIiIRKUBFRCL6/wszoLQ27BmJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(5,5))\n",
    "plt.errorbar(grid_refined, [np.mean(vals) for vals in all_scores], \n",
    "             yerr=[np.std(vals) for vals in all_scores], c=\"blue\", label=\"1. Std. Error\", alpha=0.6)\n",
    "plt.plot(grid_refined,[np.mean(vals) for vals in all_scores], 'ro', label=\"Avg. Accuracy\", markersize=9)\n",
    "# Uncomment to see standard error bounds of the best model: \n",
    "# plt.plot(grid_refined, np.repeat(best_new_accuracy+np.std(all_scores[new_chosen_k-1]), 10), \"g--\")\n",
    "# plt.plot(grid_refined, np.repeat(best_new_accuracy-np.std(all_scores[new_chosen_k-1]), 10), \"g--\")\n",
    "plt.xticks(grid_refined)\n",
    "plt.xlabel('k')\n",
    "t = plt.title('knn: Tuning k ')\n",
    "plt.ylim(0.8, 1.1)\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.legend(loc=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9333333333333333"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Performance of the best model on test set: \n",
    "knn_final = KNeighborsClassifier(n_neighbors=new_chosen_k)\n",
    "knn_final.fit(X_trainvalid, y_trainvalid)\n",
    "knn_final.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a function in scikit-learn that allows performing grid search automatically. Such function is [`GridSearchCV`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html). The inputs of the functions are: \n",
    "* estimator: the algorithm to be evaluated \n",
    "* parameter_grid: the grid over which to perform the search \n",
    "* scoring: the scoring method, as seen for the cross_val_score function\n",
    "* cv: the cross-validation option (it works exactly as in cross_val_score)\n",
    "\n",
    "After finding the best parameter set, the function works as a proper estimator; that is, we can use the `fit`, `predict` and `score` methods with it. Training and prediction will be automatically performed with the parameter set that yields the best cross-validation result. The grid must be specified as a Python dictionary, in which each key corresponds to the hyperparameter name, and each value to its grid.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9333333333333333"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "# A separate validation set is not needed anymore; we use 10% data for the test set\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.1, \n",
    "                                                    stratify=y, random_state=1)\n",
    "\n",
    "param_grid = {\"n_neighbors\" : np.arange(1, 11, 1)}\n",
    "# Perform grid search with five-fold cross validation: \n",
    "grid_search = GridSearchCV(KNeighborsClassifier(), param_grid=param_grid, cv=5)\n",
    "grid_search.fit(X_train, y_train)\n",
    "# Test set score:\n",
    "grid_search.score(X_test,y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, results match the ones of manual grid search. We can also check the best chosen parameter: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'n_neighbors': 6}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_search.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and the best cross-validation score: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9851851851851852"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_search.best_score_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to obtain the best chosen estimator (which can give us access to learned parameters, feature importances, and so on): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',\n",
       "                     metric_params=None, n_jobs=None, n_neighbors=6, p=2,\n",
       "                     weights='uniform')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_search.best_estimator_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`GridSearchCV` also gives the possibility to inspect the results within each cross-validation fold, and for each hyperparameter value. This comes in form of a dictionary, so we can turn it into a `pandas` dataframe for a better visualization: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean_fit_time</th>\n",
       "      <th>std_fit_time</th>\n",
       "      <th>mean_score_time</th>\n",
       "      <th>std_score_time</th>\n",
       "      <th>param_n_neighbors</th>\n",
       "      <th>params</th>\n",
       "      <th>split0_test_score</th>\n",
       "      <th>split1_test_score</th>\n",
       "      <th>split2_test_score</th>\n",
       "      <th>split3_test_score</th>\n",
       "      <th>split4_test_score</th>\n",
       "      <th>mean_test_score</th>\n",
       "      <th>std_test_score</th>\n",
       "      <th>rank_test_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000498</td>\n",
       "      <td>0.000156</td>\n",
       "      <td>0.001503</td>\n",
       "      <td>0.000225</td>\n",
       "      <td>1</td>\n",
       "      <td>{'n_neighbors': 1}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.000487</td>\n",
       "      <td>0.000138</td>\n",
       "      <td>0.001585</td>\n",
       "      <td>0.000216</td>\n",
       "      <td>2</td>\n",
       "      <td>{'n_neighbors': 2}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.948148</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.000578</td>\n",
       "      <td>0.000163</td>\n",
       "      <td>0.001707</td>\n",
       "      <td>0.000312</td>\n",
       "      <td>3</td>\n",
       "      <td>{'n_neighbors': 3}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.000392</td>\n",
       "      <td>0.000020</td>\n",
       "      <td>0.001440</td>\n",
       "      <td>0.000039</td>\n",
       "      <td>4</td>\n",
       "      <td>{'n_neighbors': 4}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.000410</td>\n",
       "      <td>0.000016</td>\n",
       "      <td>0.001487</td>\n",
       "      <td>0.000043</td>\n",
       "      <td>5</td>\n",
       "      <td>{'n_neighbors': 5}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>2.962963e-02</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000470</td>\n",
       "      <td>0.000094</td>\n",
       "      <td>0.001606</td>\n",
       "      <td>0.000184</td>\n",
       "      <td>6</td>\n",
       "      <td>{'n_neighbors': 6}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.985185</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.000356</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>0.001347</td>\n",
       "      <td>0.000009</td>\n",
       "      <td>7</td>\n",
       "      <td>{'n_neighbors': 7}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>2.962963e-02</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.000348</td>\n",
       "      <td>0.000014</td>\n",
       "      <td>0.001361</td>\n",
       "      <td>0.000040</td>\n",
       "      <td>8</td>\n",
       "      <td>{'n_neighbors': 8}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.970370</td>\n",
       "      <td>2.771598e-02</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.000369</td>\n",
       "      <td>0.000080</td>\n",
       "      <td>0.001296</td>\n",
       "      <td>0.000100</td>\n",
       "      <td>9</td>\n",
       "      <td>{'n_neighbors': 9}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.970370</td>\n",
       "      <td>2.771598e-02</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.000312</td>\n",
       "      <td>0.000009</td>\n",
       "      <td>0.001198</td>\n",
       "      <td>0.000038</td>\n",
       "      <td>10</td>\n",
       "      <td>{'n_neighbors': 10}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mean_fit_time  std_fit_time  mean_score_time  std_score_time  \\\n",
       "0       0.000498      0.000156         0.001503        0.000225   \n",
       "1       0.000487      0.000138         0.001585        0.000216   \n",
       "2       0.000578      0.000163         0.001707        0.000312   \n",
       "3       0.000392      0.000020         0.001440        0.000039   \n",
       "4       0.000410      0.000016         0.001487        0.000043   \n",
       "5       0.000470      0.000094         0.001606        0.000184   \n",
       "6       0.000356      0.000004         0.001347        0.000009   \n",
       "7       0.000348      0.000014         0.001361        0.000040   \n",
       "8       0.000369      0.000080         0.001296        0.000100   \n",
       "9       0.000312      0.000009         0.001198        0.000038   \n",
       "\n",
       "  param_n_neighbors               params  split0_test_score  \\\n",
       "0                 1   {'n_neighbors': 1}           0.962963   \n",
       "1                 2   {'n_neighbors': 2}           0.962963   \n",
       "2                 3   {'n_neighbors': 3}           0.962963   \n",
       "3                 4   {'n_neighbors': 4}           0.962963   \n",
       "4                 5   {'n_neighbors': 5}           0.925926   \n",
       "5                 6   {'n_neighbors': 6}           0.962963   \n",
       "6                 7   {'n_neighbors': 7}           0.925926   \n",
       "7                 8   {'n_neighbors': 8}           0.925926   \n",
       "8                 9   {'n_neighbors': 9}           0.925926   \n",
       "9                10  {'n_neighbors': 10}           0.962963   \n",
       "\n",
       "   split1_test_score  split2_test_score  split3_test_score  split4_test_score  \\\n",
       "0           0.962963           0.962963           0.962963           0.962963   \n",
       "1           0.925926           0.925926           0.962963           0.962963   \n",
       "2           0.962963           0.962963           0.962963           0.962963   \n",
       "3           0.962963           0.962963           0.962963           0.962963   \n",
       "4           1.000000           1.000000           0.962963           1.000000   \n",
       "5           1.000000           1.000000           0.962963           1.000000   \n",
       "6           1.000000           1.000000           0.962963           1.000000   \n",
       "7           0.962963           1.000000           0.962963           1.000000   \n",
       "8           0.962963           1.000000           0.962963           1.000000   \n",
       "9           0.962963           1.000000           0.962963           1.000000   \n",
       "\n",
       "   mean_test_score  std_test_score  rank_test_score  \n",
       "0         0.962963    1.110223e-16                7  \n",
       "1         0.948148    1.814437e-02               10  \n",
       "2         0.962963    1.110223e-16                7  \n",
       "3         0.962963    1.110223e-16                7  \n",
       "4         0.977778    2.962963e-02                2  \n",
       "5         0.985185    1.814437e-02                1  \n",
       "6         0.977778    2.962963e-02                2  \n",
       "7         0.970370    2.771598e-02                5  \n",
       "8         0.970370    2.771598e-02                5  \n",
       "9         0.977778    1.814437e-02                4  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "results = pd.DataFrame(grid_search.cv_results_)\n",
    "display(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each row corresponds to a different hyperparameter value; each column corresponds to various cross-validation results (estimation times, fold-specific results, average score and scores standard deviations, and final rank of the parameter value)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note on Grid Search**: if you see that the cross-validation model performance keeps improving at extreme values of the grid (for instance, the error keeps decreasing or the $f_1$ score keeps increasing), you should expand the grid with values that are beyond that extreme. \n",
    "\n",
    "In the examples of the notebook \"Cross-Validation in Practice\", we will see how to perform grid search with estimators that have more than one hyperparameter that must be tuned. You can find more on grid search in the [scikit-learn documentation](https://scikit-learn.org/stable/modules/grid_search.html#exhaustive-grid-search). \n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Random Search**. Instead of specifying a list of hyperparameter values, we can try a random search. The principle of random search is similar to the one of grid search, but in this case the algorithm searches automatically for optimal values across the hyperparameter space. \n",
    "\n",
    "In particular, we can assign a probability distribution to each hyperparameter, and let the algorithm pick values at random from such distributions. The main advantage is that you don't need to specify a grid of value for all hyperparameters, and so there is no need to explore all their possible combinations. It seems clear that this method is preferable when an algorithm has many (two or more) hyperparameters that need to be optimized. However, you need to specify the probability distributions for all hyperparameters. The parameters of such distributions, in general, will specify the main location and variability of such search. The distribution must be defined within the same range of admissible values of the hyperparameters (for example: a hyperparameter that can only be a positive real number can be assigned to a Gamma, or Exponential, distribution).\n",
    "\n",
    "In scikit-learn this can be implemented with the `RandomizedSearchCV` function. It works similarly to `GridSearchCV`, but instead of a parameter grid it takes a dictionary of probability distributions (which are available, for instance, in the `scipy.stats` module). This parameter is called `param_distributions`. Alternatively, values can be given in form of a list, in which case `RandomizedSearchCV` will draw values *without replacement*. Furthermore, the function also takes a `n_iter` argument, which specifies how many values must be drawn from the specified distribution(s). See the `scikit-learn` [documentation](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RandomizedSearchCV.html) for more information.\n",
    "\n",
    "Here, we see an example of randomized parameter search for $k$ in k-nearest-neighbors. We use 30 draws from a discrete uniform distribution over the space $k \\in \\{1,...,50\\}$ to perform the search; therefore, we specify a list of values between 1 and 50 (so that values can be drawn without replacement, as this is the only one parameter we are tuning). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9333333333333333"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "# param_grid: range of values for sampling without replacement\n",
    "param_grid = {'n_neighbors': np.arange(1, 51, 1)}\n",
    "random_search = RandomizedSearchCV(KNeighborsClassifier(), param_distributions=param_grid, \n",
    "                                   n_iter=30, cv=5, random_state=1)\n",
    "random_search.fit(X_train, y_train)\n",
    "# Test set score:\n",
    "random_search.score(X_test,y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function returns the same types of values as `GridSearchCV`: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'n_neighbors': 5}"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random_search.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize the values drawn for $k$: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3,  4,  5, 11, 14, 15, 18, 20, 22, 23, 24, 25, 27, 28, 30, 31, 32,\n",
       "       33, 34, 36, 37, 39, 40, 41, 42, 43, 46, 47, 49, 50])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = pd.DataFrame(random_search.cv_results_)\n",
    "np.sort(np.unique([list((results[\"params\"][x].values()))[0] for x in range(30)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the results: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean_fit_time</th>\n",
       "      <th>std_fit_time</th>\n",
       "      <th>mean_score_time</th>\n",
       "      <th>std_score_time</th>\n",
       "      <th>param_n_neighbors</th>\n",
       "      <th>params</th>\n",
       "      <th>split0_test_score</th>\n",
       "      <th>split1_test_score</th>\n",
       "      <th>split2_test_score</th>\n",
       "      <th>split3_test_score</th>\n",
       "      <th>split4_test_score</th>\n",
       "      <th>mean_test_score</th>\n",
       "      <th>std_test_score</th>\n",
       "      <th>rank_test_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.000581</td>\n",
       "      <td>0.000201</td>\n",
       "      <td>0.001754</td>\n",
       "      <td>0.000135</td>\n",
       "      <td>3</td>\n",
       "      <td>{'n_neighbors': 3}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000397</td>\n",
       "      <td>0.000030</td>\n",
       "      <td>0.001448</td>\n",
       "      <td>0.000082</td>\n",
       "      <td>4</td>\n",
       "      <td>{'n_neighbors': 4}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>0.000305</td>\n",
       "      <td>0.000005</td>\n",
       "      <td>0.001153</td>\n",
       "      <td>0.000028</td>\n",
       "      <td>5</td>\n",
       "      <td>{'n_neighbors': 5}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>2.962963e-02</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>0.000299</td>\n",
       "      <td>0.000001</td>\n",
       "      <td>0.001158</td>\n",
       "      <td>0.000008</td>\n",
       "      <td>11</td>\n",
       "      <td>{'n_neighbors': 11}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>0.000313</td>\n",
       "      <td>0.000024</td>\n",
       "      <td>0.001182</td>\n",
       "      <td>0.000016</td>\n",
       "      <td>14</td>\n",
       "      <td>{'n_neighbors': 14}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.970370</td>\n",
       "      <td>1.481481e-02</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mean_fit_time  std_fit_time  mean_score_time  std_score_time  \\\n",
       "4        0.000581      0.000201         0.001754        0.000135   \n",
       "5        0.000397      0.000030         0.001448        0.000082   \n",
       "25       0.000305      0.000005         0.001153        0.000028   \n",
       "29       0.000299      0.000001         0.001158        0.000008   \n",
       "19       0.000313      0.000024         0.001182        0.000016   \n",
       "\n",
       "   param_n_neighbors               params  split0_test_score  \\\n",
       "4                  3   {'n_neighbors': 3}           0.962963   \n",
       "5                  4   {'n_neighbors': 4}           0.962963   \n",
       "25                 5   {'n_neighbors': 5}           0.925926   \n",
       "29                11  {'n_neighbors': 11}           0.962963   \n",
       "19                14  {'n_neighbors': 14}           0.962963   \n",
       "\n",
       "    split1_test_score  split2_test_score  split3_test_score  \\\n",
       "4            0.962963           0.962963           0.962963   \n",
       "5            0.962963           0.962963           0.962963   \n",
       "25           1.000000           1.000000           0.962963   \n",
       "29           0.962963           1.000000           0.962963   \n",
       "19           0.962963           1.000000           0.962963   \n",
       "\n",
       "    split4_test_score  mean_test_score  std_test_score  rank_test_score  \n",
       "4            0.962963         0.962963    1.110223e-16                6  \n",
       "5            0.962963         0.962963    1.110223e-16                6  \n",
       "25           1.000000         0.977778    2.962963e-02                1  \n",
       "29           1.000000         0.977778    1.814437e-02                2  \n",
       "19           0.962963         0.970370    1.481481e-02                4  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(results.sort_values([\"param_n_neighbors\"]).head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For completeness, it is shown here how to perform random search with a probability distribution. Here, we use the discrete uniform distribution provided with the `scipy.stats.randint` [function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.randint.html). 20 draws will be performed.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9333333333333333"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from  scipy import stats\n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "param_grid = {'n_neighbors': stats.randint(1, 50)}\n",
    "random_search = RandomizedSearchCV(KNeighborsClassifier(), param_distributions=param_grid, \n",
    "                                   n_iter=20, cv=5, random_state=1)\n",
    "random_search.fit(X_train, y_train)\n",
    "# Test set score:\n",
    "random_search.score(X_test,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean_fit_time</th>\n",
       "      <th>std_fit_time</th>\n",
       "      <th>mean_score_time</th>\n",
       "      <th>std_score_time</th>\n",
       "      <th>param_n_neighbors</th>\n",
       "      <th>params</th>\n",
       "      <th>split0_test_score</th>\n",
       "      <th>split1_test_score</th>\n",
       "      <th>split2_test_score</th>\n",
       "      <th>split3_test_score</th>\n",
       "      <th>split4_test_score</th>\n",
       "      <th>mean_test_score</th>\n",
       "      <th>std_test_score</th>\n",
       "      <th>rank_test_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>0.000305</td>\n",
       "      <td>4.125675e-06</td>\n",
       "      <td>0.001094</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>1</td>\n",
       "      <td>{'n_neighbors': 1}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.110223e-16</td>\n",
       "      <td>11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0.000300</td>\n",
       "      <td>4.101908e-07</td>\n",
       "      <td>0.001167</td>\n",
       "      <td>0.000092</td>\n",
       "      <td>2</td>\n",
       "      <td>{'n_neighbors': 2}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.948148</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.000329</td>\n",
       "      <td>1.130718e-05</td>\n",
       "      <td>0.001260</td>\n",
       "      <td>0.000051</td>\n",
       "      <td>6</td>\n",
       "      <td>{'n_neighbors': 6}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.985185</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>0.000306</td>\n",
       "      <td>9.250641e-06</td>\n",
       "      <td>0.001251</td>\n",
       "      <td>0.000140</td>\n",
       "      <td>7</td>\n",
       "      <td>{'n_neighbors': 7}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>2.962963e-02</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>0.000303</td>\n",
       "      <td>6.021379e-06</td>\n",
       "      <td>0.001223</td>\n",
       "      <td>0.000165</td>\n",
       "      <td>8</td>\n",
       "      <td>{'n_neighbors': 8}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.970370</td>\n",
       "      <td>2.771598e-02</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.000380</td>\n",
       "      <td>2.523428e-05</td>\n",
       "      <td>0.001450</td>\n",
       "      <td>0.000131</td>\n",
       "      <td>9</td>\n",
       "      <td>{'n_neighbors': 9}</td>\n",
       "      <td>0.925926</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.970370</td>\n",
       "      <td>2.771598e-02</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.000363</td>\n",
       "      <td>1.286754e-05</td>\n",
       "      <td>0.001371</td>\n",
       "      <td>0.000021</td>\n",
       "      <td>10</td>\n",
       "      <td>{'n_neighbors': 10}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.000354</td>\n",
       "      <td>8.991131e-06</td>\n",
       "      <td>0.001338</td>\n",
       "      <td>0.000033</td>\n",
       "      <td>12</td>\n",
       "      <td>{'n_neighbors': 12}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.000482</td>\n",
       "      <td>6.869080e-05</td>\n",
       "      <td>0.001688</td>\n",
       "      <td>0.000081</td>\n",
       "      <td>13</td>\n",
       "      <td>{'n_neighbors': 13}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0.000370</td>\n",
       "      <td>7.939970e-05</td>\n",
       "      <td>0.001219</td>\n",
       "      <td>0.000067</td>\n",
       "      <td>13</td>\n",
       "      <td>{'n_neighbors': 13}</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.962963</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.977778</td>\n",
       "      <td>1.814437e-02</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    mean_fit_time  std_fit_time  mean_score_time  std_score_time  \\\n",
       "8        0.000305  4.125675e-06         0.001094        0.000004   \n",
       "10       0.000300  4.101908e-07         0.001167        0.000092   \n",
       "6        0.000329  1.130718e-05         0.001260        0.000051   \n",
       "14       0.000306  9.250641e-06         0.001251        0.000140   \n",
       "12       0.000303  6.021379e-06         0.001223        0.000165   \n",
       "3        0.000380  2.523428e-05         0.001450        0.000131   \n",
       "4        0.000363  1.286754e-05         0.001371        0.000021   \n",
       "5        0.000354  8.991131e-06         0.001338        0.000033   \n",
       "2        0.000482  6.869080e-05         0.001688        0.000081   \n",
       "11       0.000370  7.939970e-05         0.001219        0.000067   \n",
       "\n",
       "   param_n_neighbors               params  split0_test_score  \\\n",
       "8                  1   {'n_neighbors': 1}           0.962963   \n",
       "10                 2   {'n_neighbors': 2}           0.962963   \n",
       "6                  6   {'n_neighbors': 6}           0.962963   \n",
       "14                 7   {'n_neighbors': 7}           0.925926   \n",
       "12                 8   {'n_neighbors': 8}           0.925926   \n",
       "3                  9   {'n_neighbors': 9}           0.925926   \n",
       "4                 10  {'n_neighbors': 10}           0.962963   \n",
       "5                 12  {'n_neighbors': 12}           0.962963   \n",
       "2                 13  {'n_neighbors': 13}           0.962963   \n",
       "11                13  {'n_neighbors': 13}           0.962963   \n",
       "\n",
       "    split1_test_score  split2_test_score  split3_test_score  \\\n",
       "8            0.962963           0.962963           0.962963   \n",
       "10           0.925926           0.925926           0.962963   \n",
       "6            1.000000           1.000000           0.962963   \n",
       "14           1.000000           1.000000           0.962963   \n",
       "12           0.962963           1.000000           0.962963   \n",
       "3            0.962963           1.000000           0.962963   \n",
       "4            0.962963           1.000000           0.962963   \n",
       "5            0.962963           1.000000           0.962963   \n",
       "2            0.962963           1.000000           0.962963   \n",
       "11           0.962963           1.000000           0.962963   \n",
       "\n",
       "    split4_test_score  mean_test_score  std_test_score  rank_test_score  \n",
       "8            0.962963         0.962963    1.110223e-16               11  \n",
       "10           0.962963         0.948148    1.814437e-02               13  \n",
       "6            1.000000         0.985185    1.814437e-02                1  \n",
       "14           1.000000         0.977778    2.962963e-02                2  \n",
       "12           1.000000         0.970370    2.771598e-02                8  \n",
       "3            1.000000         0.970370    2.771598e-02                8  \n",
       "4            1.000000         0.977778    1.814437e-02                3  \n",
       "5            1.000000         0.977778    1.814437e-02                3  \n",
       "2            1.000000         0.977778    1.814437e-02                3  \n",
       "11           1.000000         0.977778    1.814437e-02                3  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Display the results of RandomizedSearchCV(); as you can see, some parameter values are duplicated \n",
    "results = pd.DataFrame(random_search.cv_results_)\n",
    "display(results.sort_values([\"param_n_neighbors\"]).head(10))"
   ]
  }
 ],
 "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}