{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LezvU_JFKPIb"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sklearn.linear_model as linear_model\n",
    "from sklearn.model_selection import train_test_split\n",
    "import sklearn.preprocessing as preprocessing\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AIOhMmjKKTPs"
   },
   "source": [
    "# 1 Polynomial Regression with SGD\n",
    "\n",
    "In this exercise we will have a look again at overfitting and underfitting using the example of a polynomial regression.\n",
    "However, we will find the optimal coefficient not by solving the linear problem as in Exercise 2, but by minimizing a loss function through Stochastic Gradient Descent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function generates data from a polynomial of degree $d$ with given coefficients $\\{c_k\\}_{k=0}^{d}$ and noise level from a Gaussian of standard deviation $\\xi$, namely $y_i = \\sum_{k=0}^{d} c_k x_i^k + \\xi \\mathcal{N}(0, 1)$ where $x_i$ are uniformly distributed points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "_xLwEsVtKTEL"
   },
   "outputs": [],
   "source": [
    "# This function generates data from a polynome, possibly with noise\n",
    "def generate_polynome_data(polynome, noise_std, n, bound, seed):\n",
    "    \"\"\"\n",
    "    arguments:\n",
    "        - polynome  : list of coefficients representing the polynom (polynome[i] multiplies X^i)\n",
    "        - noise_std : standard deviation of the gaussian noise\n",
    "        - n         : number of data points\n",
    "        - bound     : for points x the sampling inverval is +-bound\n",
    "    returns:\n",
    "        - (X, y)    : list of n points with y = poly(x) + gaussian noise\n",
    "    \"\"\"\n",
    "    ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fYggcS-iUg1q"
   },
   "source": [
    "## 1.1 Generate data\n",
    "\n",
    "Generate $n$ noisy data points (say $5 \\leqslant n \\leqslant 10$) from the polynomial $y_i = x_i^2  + 0.5 x_i^3$ with a given noise level.\n",
    "These represent the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "MrL0B_lUKTBq"
   },
   "outputs": [],
   "source": [
    "# Create 100 noisy datapoints within the training interval [-3, 3]\n",
    "polynom = np.array([0, 0, 1, 0.5])  # function y = x^2 + 0.5 * x^3\n",
    "X_train, y_train = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 279
    },
    "id": "INC5GWlNQvAT",
    "outputId": "208d076c-e89c-4331-b333-e745599519dd"
   },
   "outputs": [],
   "source": [
    "# Show the training data in a plot\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Yaaz_XX1bZi-"
   },
   "source": [
    "We want to create some data for the test set too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NFTRqawrO30F"
   },
   "outputs": [],
   "source": [
    "X_test, y_test = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zBQF-qFHU3QB"
   },
   "source": [
    "## 1.2 Polynomials of different degrees\n",
    "The loss function that we use for the training is the L2 metrics\n",
    "$L(c) = \\frac{1}{2n}\\sum_{x_i \\in X}(y_i-\\hat{y_i})^2$, where $\\hat{y_i}$ is the output of a polynomial on the input $x_i$ and $y_i$ is the ground-truth output. \n",
    "Note: add a L2 regularization term to the loss function. \n",
    "\n",
    "When you generate the features, rescale the powers with the factorial of the degree, namely generate $x_i^k/k!$ instead of $x_i^k$ as done in the Exercise 2. \n",
    "This would avoid overflow in the fitting procedure. \n",
    "\n",
    "What are the gradients with respect to the coefficients $c$?\n",
    "\n",
    "$\\frac{\\partial L(c)}{\\partial c_k} = - \\frac{1}{n}\\sum_{x_i \\in X} (y_i-\\hat{y_i}) x_i^k$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to write the functions to compute the loss function and its gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function that generates the feature matrix\n",
    "def create_features(X, d):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to predict the labels\n",
    "def predict_label(X, c):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to compute the loss\n",
    "def loss(X, y, c, lambd):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CYACZ6d_Y2Nr"
   },
   "outputs": [],
   "source": [
    "# Fucntion to compute the gradient of the loss\n",
    "def grad(X, y, c, lambd):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write a function that performs the training using the Stochastic Gradient Descent scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 370
    },
    "id": "dmt7_oasOLMN",
    "outputId": "50e73e96-303e-4d7b-cbb4-e163695e995d"
   },
   "outputs": [],
   "source": [
    "# Perform the SGD training\n",
    "def sgd_training(\n",
    "    X_train,\n",
    "    y_train,\n",
    "    X_test,\n",
    "    y_test,\n",
    "    d,\n",
    "    BATCHSIZE=16,\n",
    "    EPOCHS=1000,\n",
    "    lr=0.001,\n",
    "    lambd=1e-12,\n",
    "):\n",
    "    # create the features\n",
    "    ...\n",
    "\n",
    "    # set initial values for the parameters\n",
    "    ...\n",
    "\n",
    "    # for several epochs (= runs over complete training data)\n",
    "    for e in range(EPOCHS):\n",
    "        # compute current errors on test and train\n",
    "        ...\n",
    "\n",
    "        # split the training set into batches and do a gradient step for each batch\n",
    "        ...\n",
    "\n",
    "        # run a gradient step on every batch\n",
    "        ...\n",
    "\n",
    "    return ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the training and test errors for several degree of polynomials\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the learning curves\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us compute the final test error as a function of the model complexity (degree of the polynomial). \n",
    "Show that you recover the well-known U-shaped curve of the bias-variance tradeoff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the best MSE on the test set\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot the curves of the different polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the learned polynomials\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2:  SGD for Linear Regression on a health dataset\n",
    "\n",
    "In this exercise we want to repeat how to load a real-world dataset.\n",
    "\n",
    "Then, we want to implement SGD to infer the parameters and compare this with the LinearRegression classes results from sklearn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import diabetes dataset from scikit-learn toy-dataset examples\n",
    "from sklearn.datasets import load_diabetes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load boston dataset\n",
    "diabetes = load_diabetes()\n",
    "\n",
    "# Check what it contains\n",
    "print(diabetes.DESCR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the dimension of X is (n, d)= (442, 10) the dimension of y is (n,)= (442,)\n"
     ]
    }
   ],
   "source": [
    "# Extract the dataset\n",
    "X = diabetes.data\n",
    "y = diabetes.target\n",
    "# print the dimension of the X matrix and the y vector\n",
    "print(\"the dimension of X is (n, d)=\", X.shape, \"the dimension of y is (n,)=\", y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAALICAYAAABiqwZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZRlVX328e/Tg7SIgAoqTrTiAAalA+WAosEpGmMUFUTAATW2s68ah86r0U7emNVGE2eNrSIYfQkqigRU8EUZ7ATpAoEGFDGCUYIoo0wNTffv/eOelm1RXbequ6tuDd/PWnfVuefss/fvnHaxnrvd99xUFZIkSZJ65g26AEmSJGk6MSBLkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJDQOyJEmS1Fgwkcbrrvr5tPrZvYU7PSSDrkGSJEmzy4QCMhvWT1IZkiRJ0vQwsYC8/vZJKkMzydDQ0DJg0aDr0J2sHR4eXjHoIiRJmukmFJDLgKyeRcPDw8sHXYT+0NDQ0PJB1yBJ0mwwwRnkdZNUhiRJkjQ9uMRCkiRJarjEQpIkSWq4xEKSJElquMRCkiRJahiQJUmSpMbE1iBvcImFJEmSZjdnkCVJkqSGX9KTJEmSGs4gS5IkSY2JBeTbDciSJEma3Sb4QyEusZAkSdLs5hILSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGrP6OchJngV8FJgPfK6qVow4fjjwQeDybtcnqupzU1qkJEmSppUJziCvn6Qytr4k84FPAs8AfgWsTnJ8VV00oukxVfXGKS9QkiRJ09JsXmLxWOBnVfVzgCT/BjwPGBmQJUmSpN+b0TPISZYCS5tdK6tqZbd9f+CXzbFfAY8bpZsXJnky8FPgrVX1y1HaSJIkaY6Y0WuQuzC8sm/DTft34OiqujXJa4CjgKduleIkSZI0I00oINc0m0Hu43Lggc37B3DHl/EAqKqrm7efA/5xCuqSJEnSNDab1yCvBh6W5MH0gvGLgUPbBkl2qaorurfPBX48tSVKkiRpupngEouZM4NcVbcneSNwEr3HvB1RVRcm+TtguKqOB96c5LnA7cA1wOEDK1iSJEnTwoz+kl4/VfUt4Fsj9r232f5r4K+nui5JkiRNX7N2BlmSJEnaHLN6BlmSJEmaqIk9xcIZZEmSJM1yziBLkiRJDdcgS5IkSQ1nkCVJkqSGa5AlSZKkhjPIkiRJUmOCM8gbJqsOSZIkaVrwS3qSJElSY2IBeUNNUhmSJEnS9DCxJRbrXWIhSZKk2W2CSywMyJIkSZrd/JKeJEmS1JhgQHYNsiRJkma3CS6xMCBLkiRpdnMGWZIkSWoYkCVJkqTGBAPyZJUhSZIkTQ8GZEmSJKkxoYC8wYAsSZKkWW6Cv6SXyapDkiRJmhYmOINsQNbMMDQ0tAxYNOg6ptjioaGh5YMuYgDWDg8Prxh0EZKk2WNiAdkZZM0ci4aHh5cPughNvjn6oUCSNImcQZYkSZIaziBLkiRJDQOyJEmS1JhgQJ43WXVIkiRJ08KEAvJ6Z5ClaWuOPrkD5ubTO3xyhyRNImeQpdljTj65Y45+MNh9Dn4oAD8YSJoirkGWNNPNyQ8Gc9Ec/VAgaQAmtsRigzPIkiRJmt0mGJCdQZYkSdLsNqEp4Q0bMq1e/SR5VpKLk/wsybJRjm+T5Jju+A+TLJ7I/ZAkSdLsM2uXWCSZD3wSeAbwK2B1kuOr6qKm2auAa6vqoUleDHwAOHjqq5UkSdJ0MbGAXDNqicVjgZ9V1c8Bkvwb8DygDcjPA5Z3218DPpEkVVVTWagkSZKmjxk9g5xkKbC02bWyqlZ22/cHftkc+xXwuBFd/L5NVd2e5HrgXsBVk1OxJEmSpruJBWSm1wxyF4ZX9m0oSZIkjdOEAvLtM2uJxeXAA5v3D+j2jdbmV0kWADsAV09NeZIkSZqOZvQMch+rgYcleTC9IPxi4NARbY4HXg78J3Ag8D3XH0ua7uborweCPysuaYrM2oDcrSl+I3ASMB84oqouTPJ3wHBVHQ98HvjXJD8DrqEXoiVpuvPXA+eIOfiBQJoWZm1ABqiqbwHfGrHvvc32WuCgqa5LkiRJ09fE1iBnZgVkSZI2ZYYsVZkJy0pcBqJZZ4IzyJIkzRpzbqnKJH0o2H0rh3gDtwZuYgHZGWRJkmayaf+hYAbMmGsOmOASi8kqQ5IkSZoeZvWX9CRJkqSJcgZZkiRJakxwDfJklSFJs9ckfDFqMp5s4BejJKkzsRnkyapCkmY3vxglSTOIM8iSJElSY0IBecNkVSFJkiRNE84gS5IkSQ3XIEuSJEkNZ5AlSTOCTwORNFUm+EMhkiQNjE8DkTQlJrjEoiarDkmSJGlacAZZkiRJakzwp6adQZYkSdLs5gyyJEmS1JhgQHYGWZIkSbObX9KTJEmSGi6xkCRJkhousZAGxB89kCRpejIgS4Pjjx5IkjQNuQZZkiRJajiDLEmSJDUMyJIkSVJjYgG5DMiSJEma3SY4g7xhsuqQJEmSpgW/pCdJkiQ1XIMsSZIkNSa4BtklFpIkSZrdnEGWJEmSGs4gS5IkSY05O4Oc5J7AMcBi4DLgRVV17Sjt1gNrurf/XVXPnaoaJUmSNPXm8gzyMuCUqlqRZFn3/l2jtLulqpZMbWmSJEkalLn8HOTnAft320cBpzJ6QJYkSdIckprAr+M96f5Pm1ZrLM64/JRs7rlJrquqHbvtANdufD+i3e3AucDtwIqqOm5zxxyEdVf9fFr9m43m5fv81aBL6OudrB90CX2dUtsPuoS+drtt+n/IPvCa0wZdQl/P32Vo0CX0tYj5gy6hr9tmwKTPfhu2G3QJY/rV/On/30aAD1x29GbnBc1NE/yhkOn1H5MkS4Glza6VVbWyOf7/gPuOcuq72zdVVUk2FSR3rarLkzwE+F6SNVX1X1tauyRJkqanGb0GuQvDK8c4/vRNHUtyZZJdquqKJLsAv9lEH5d3f3+e5FTgjwEDsiRJ0iw1byKNN9SGafXaQscDL++2Xw58c2SDJPdIsk23vRPwROCiLR1YkiRJ09eMnkHeQiuAryR5FfAL4EUASYaA11bVXwJ7AJ9JsoHeh4kVVWVAliRJmsXm7FMsqupq4Gmj7B8G/rLb/g/gUVNcmiRJkgZoLs8gS5IkSXdiQJYkSZIaBmRJkiSpYUCWJEmSGhMMyDPjF3MkSZKkzeUMsiRJktSYWEDeYECWJEnS7DahgFzUZNUhSZIkTQsusZAkSZIaLrGQJEmSGs4gS5IkSQ1nkCVJkqSGM8iSJElSY95EGm+oDdPqJUmStKWSLE5ywSj7/y7J0/ucuzzJ2zdz3M0+d4w+v5Vkx+71+s3sY/8kJ3Tbz02ybGvWOIE6Fic5dBzt7pfka1tz7AnNIG9wiYUkSZojquq9g65hoqrq2dALl8DrgU9tYX/HA8dvcWGbkCRAqkad+VwMHAr837H6qKr/AQ7cmnVNKCCvu+3yAAwNDS0DFm3NQjbH0NDQ8kHXMN3953e+MugSJEkaqCRLgaXNrpVVtXJEs/lJPgs8AbgceB7waeCEqvpakmcD/wzcBKwCHlJVz+nOfWSSU4EHAR+pqo+NUcu7gZcDvwF+CZzd7d8N+CSwM3Az8Oqq+kmSI4HfAUPAfYF3dvXsAhwDbE8vz72uqs5IclnXdgWwW5Jzge8C9wG+XlXHdeN9GfhKVX2zz707HBiqqjd2tazt+t8eeFtVnZBkUXevhoDbu/3f7859PrADcH/gS1X1t114Pwn4IbAP8OwkbwT+DCjg76vqmO4a9uiu4SjgY92+/YFtgE9W1We6/k6oqj27MZ8LbAvsBnyjqt6ZZD7w+a7GAo6oqg9v6ronFJAbi4aHh5dv5rmaQuuu+vn7Bl2DJEmD1IXhkYF4pIcBh1TVq5N8BXjhxgNdAPwM8OSqujTJ0SPO3R14CnB34OIkn66qdSMHSLIP8GJgCb0Mdg5dQO7qe21VXZLkcfRmfp/aHdsF2K8b53jga/RmVk+qqvd34W/bEcMtA/asqiXd2H8CvBU4LskO9D4IvLzPPRnNYuCx9MLn95M8FHgDUFX1qCS7AycneXjX/rHAnvRC/+okJwJX0bvfL6+qM5O8sLsnewE7de1O767h7Rs/iHQfdK6vqsck2QZYleRkuNMv2S0B/hi4ld6/x8eBewP3r6o9u752HOsiJ7QGWZIkaZa6tKrO7bbPphcEN9od+HlVXdq9HxmQT6yqW6vqKnozw/fZxBhPojejeXNV/Y5u6UKS7egF1q92s6WfoReKNzquqjZU1UVN36uBVyRZDjyqqm4Y6+Kq6jTgYUl2Bg4Bjq2q28c6ZxO+0tVyCfBzevdmP+BL3Tg/AX4BbAzI362qq6vqFuDrXVuAX1TVmd32fsDRVbW+qq4ETgMeM8rYfwq8rLtHPwTuRS9oj3RKVV1fVWuBi4Bdu1ofkuTjSZ5Fb1Z+kwzIkiRJvdnGjdYzsf+XfUvOhV4eu66qljSvPTbRfwCq6nTgyfSWgxyZ5GXjGOeLwEuAVwBHTLDGjUbO1o58P972N23G2AHe1NyjB1fVyaO0u9O/R1VdS2+G+lTgtcDnxhrIgCxJkjS2i+nNPi7u3h+8mf2cDhyQ5K5J7g78BUA3m3xpkoOg98W1JHuN1VGSXYErq+qz9MLe3iOa3EBvyUfrSOAt3ZgXbeY1HJRkXrdm+iH07s0ZwGFdXQ+ntxb74q79M5LcM8ldgQPord8e6Qzg4CTzuxnuJwNnjXINJwGvS7Jw41hJ7jaeopPsBMyrqmOB93Dn+/UHNncNsiRJ0pxQVbd0j0z7TpKb6C1v2Jx+zklyDHAevaUYbT+HAZ9O8h5gIfBvXbtN2R94R5J1wI3AH8wgV9XVSVZ1j6/7dlW9o6quTPJj4LjNqb/z3/TC6/b01kyvTfKprvY19L6kd3hV3dp7QAVnAccCD6D3Jb3h5oPGRt8A9u2ut+h9EfHXSa4G1ic5j164/yi9pS/ndE+/+C290D0e9we+kGTj5PBfj9U4Vf1mxu9saGhouV/SmxnWXfXzif8DT7GX7/NXgy6hr3eyftAl9HVKbT/oEvra7bbp/6jIA685bdAl9PX8XYYGXUJfi5g/6BL6uo3p/7/H/TZsN+gSxvSr+dP/v40AH7js6GxpH0m2q6obu2D2SeCSsZ6CMB0l2RZYA+xdVddvxvlH0j3VY5ztD6d7AsZExxo0l1hIkiT19+ruy2EX0nts2WcGXM+EdD948mPg45sTjucal1hIkiT10c0Wj2vGOMm9gFNGOfS0qrp6qxY2TlX1/+g9zeH3kjwT+MCIppdW1fM30cfhExzzSHpLI2YcA7IkSdJW1IXgJYOuo5+qOoneF980gkssJEmSpIYBWZIkSWoYkCVJkqSGAVmSJElqGJAlSZKkhgFZkiRJahiQJUmSpIYBWZIkSWoYkCVJkqSGAVmSJElq+FPTs9zL9/mrQZfQ11Fn/9OgS+jrvCVvG3QJfR103ysGXUJfCxdtGHQJfV07tPegS+jrulU3DrqEvm66fptBl9DXLWsXDrqEvnYd+uWgSxjThpsHXYE0OZxBliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSNKcluVuSE5Ocl+SCJAcn2SfJaUnOTnJSkl2S7JDk4iSP6M47OsmrB12/tr4Fgy5gkIaGhpYBiwZdx2R6OA8cdAmSJA1UkqXA0mbXyqpa2bx/FvA/VfXnXfsdgG8Dz6uq3yY5GHh/Vb0yyRuBI5N8FLhHVX12ii5DU2hOB2Rg0fDw8PJBFzGZDt31+e8bdA2SJA1SF4ZXjtFkDfBPST4AnABcC+wJfDcJwHzgiq6v7yY5CPgksNdk1q3BmesBWZIkzXFV9dMkewPPBv4e+B5wYVXtO7JtknnAHsDNwD2AX01lrZoarkGWJElzWpL7ATdX1ZeADwKPA3ZOsm93fGGSP+qavxX4MXAo8IUkCwdRsyaXM8iSJGmuexTwwSQbgHXA64DbgY9165EXAB9Jcjvwl8Bjq+qGJKcD7wFczjjLGJAlSdKcVlUnASeNcujJo+zboznvbZNWlAbKJRaSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNRYMugBNrneyftAl9HXekrcNuoS+9jr3nwddwqxw05teNegS+lqw3+MGXUJfOw3dMugS+trum6sGXUJf62+sQZfQV+Zl0CWMaZtHbDfoEqRJ4QyyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkqQ5K8mOSV7fp83iJIeOo6/FSS7YetUNRpLDk3xigudclmSnzRzrfhM9b4z+Tk0ytKX9GJAlSdJctiMwZkAGFgN9A/JkSjJ/kONPosOBCQXkJAsmp5Q7GJAlSdJctgLYLcm5ST7YvS5IsibJwU2bJ3Vt3trNFJ+R5Jzu9YTxDJTkj5Kc1fVzfpKHdftf0uz/zMYwnOTGJP+U5Dxg33aWNslQklO77eVJjupq+kWSFyT5x+4avpNk4Rg1PSbJfyQ5r6vh7t2h+3XnXpLkH5v2h3T9XpDkA5vo807X072ObO7tW5McCAwBX+7a3jXJPklOS3J2kpOS7NL1eWqSjyQZBv5Xkqcl+VHX1xFJthnPv8F4TXoC15YbGhpaBizanHM/N7EPZZIkzTpJlgJLm10rq2plt70M2LOqliR5IfBaYC9gJ2B1ktO7Nm+vqud0/W0LPKOq1nYh92h6Qa+f1wIfraovJ7kLMD/JHsDBwBOral2STwGHAV8E7gb8sKr+qht3rL53A54CPBL4T+CFVfXOJN8A/hw4bpT7chfgGODgqlqdZHvglu7wEuCPgVuBi5N8HFgPfADYB7gWODnJAVV1XNPnpq7nQuD+VbVn127HqrouyRu7ezvcBfmPA8+rqt92H1DeD7yy6/4uVTWUZBFwCfC0qvppki8CrwM+0uf+j5sBeWZYNDw8vHxzTjx31+e+byvXIknSjNKF4ZV9G8J+wNFVtR64MslpwGOA341otxD4RJIl9ELjw8dZyn8C707yAODrVXVJkqfRC5yruwB8V+A3Xfv1wLHj7PvbXSBdA8wHvtPtX0NvichoHgFcUVWrAarqd/D7IH5KVV3fvb8I2BW4F3BqVf222/9l4Mn8Yfje1PX8O/CQLmifCJy8iXr2BL7bnTsfuKI5fkzT7tKq+mn3/ijgDRiQJUmSBuatwJX0ZprnAWvHc1JV/d8kP6Q3o/utJK8BAhxVVX89yilru7C+0e3csTx25P+zfGs3xoYk66qquv0b2Ly8d2uzvX4CfWzyepLsBTyT3kz6i7hjZrg998Kq2ncTfd80zhq2mGuQJUnSXHYDsHHd7RnAwd162Z3pzY6eNaINwA70Zl43AC+lN9PZV5KHAD+vqo8B3wQeDZwCHJjk3l2beybZdRNdXEZvdhbgheO7vDFdDOyS5DHd2Hfv8wW4s4A/SbJTt076EOC0EW1GvZ5u7fS8qjoWeA+wd9e+vbcXAzsn2bc7d2GSP9pE3YuTPLR7/9JR6tgiziBLkqQ5q6quTrIqvcezfRs4HzgPKOCdVfXrJFcD67svyx0JfAo4NsnL6C1lGO/M5ouAlyZZB/wa+IequibJe+it550HrKO3XOAXo5z/t8Dnk/wf4NTNu+I7VNVt3Trfjye5K731x08fo/0VSZYB36c323tiVX1zRJuLNnE9twBf6PYBbJxhPhL4lyS3APsCBwIfS7IDvZz6EXrrl9sx1iZ5BfDVLtCvBv5lc+/DaAzIkiRpTquqkY9we8eI4+uAp45o8+hm+11du8voraHd1Dgr6D0RY+T+Y7hjfW27f7sR789glPXOVbV8U+eNPDbKuauBx4/YfWT32tjmOc320fS+lDiyn8XN9qjXwx2zxu15x/KH66zPpTdzP7Ld/iPen0LvS4RjtttcLrGQJEmSGs4gS5IkbUVJnknvcWitS6vq+YOoB6B73NuDR+x+V1WdNIh6pjsDsiRJ0lbUhc5pFTwHGc5nIpdYSJIkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJjQWbed7aoaGh5VuzkAFZvKkDQ0NDy4BFU1fKmBZv7omn1PZbsYzJcdB9rxh0CZoi63+3YdAl9JUHPWLQJfS17sgjBl1CX1mQQZfQ17obpn+N62+pQZcwpvk73jroEqRJsVkBeXh4eMXWLmQQ+oT8RcPDw2MdnzKz5MOIJEnSjOASC0mSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZI0pyW5cQrGWJ7k7Vu5z28l2XGM429Jsu1m9n1j9/d+Sb62uTVuqST/e5ztxrwXE2VAliRJmoGq6tlVdd0YTd4CbFZAbsb4n6o6cEv66CfJ/DEOjysgj+NeTIgBWZIkCUjPB5NckGRNkoO7/Z9M8txu+xtJjui2X5nk/WP09+4kP03yA+ARzf7dknwnydlJzkiye7f/yCQfS/IfSX6e5MBu/y5JTk9yblfbk7r9lyXZKcndkpyY5Lzu+MFJ3gzcD/h+ku93tX6kqeHVST48jnuyOMkF3fbhSb6Z5NQklyR5X9Pubd3YFyR5S3PuT5J8OcmPk3xt44x2V/sHkpwDHJTkkO6eX5DkA12bFcBdu+v+crfvJUnO6vZ9ZmO4bu7F4m6szya5MMnJSe7atXlzkouSnJ/k38a67gX9boymhbVDQ0PLN+fEQ9h9K5ciSdLMkmQpsLTZtbKqVo7S9AXAEmAvYCdgdZLTgTOAJwHHA/cHdunaPwkYNWgl2Qd4cdffAuAc4OyN4wOvrapLkjwO+BTw1O7YLsB+wO7deF8DDgVOqqr3d4Fw5Kzws4D/qao/78beoaquT/I24ClVdVWS7YB3J3lHVa0DXgG8ZpM3bdMeC+wJ3NzdnxOB6vp7HBDgh0lOA66l98HgVVW1qvtg8XrgQ11fV1fV3knuB5wJ7NOdc3KSA6pqWZI3VtWS7rr2AA4GnlhV65J8CjgM+OKIGh8GHFJVr07yFeCFwJeAZcCDq+rWfssxDMgzwPDw8IrNPfefHvSS9/VvJUnS7NWF4dEC8Uj7AUdX1Xrgyi7kPYZeQH5LkkcCFwH3SLILsC/w5k309STgG1V1M0CS47u/2wFPAL6aZGPbbZrzjquqDcBFSe7T7VsNHJFkYXf83BFjrQH+qZt5PaGqzhjlHtyY5HvAc5L8GFhYVWvGcU9G+m5VXd1dy9fp3bPqrvWmZv/GDxS/rKpV3blfone/NgbkY7q/jwFOrarfdud/GXgycNyIsZ9GL0Sv7u7dXYHfjFLjpc09OhtY3G2fD3w5yXGj9P0HXGIhSZI0hqq6HNiR3kztxhnlFwE3VtUNE+xuHnBdVS1pXns0x29tttONfzq9wHg5cGSSl42o76fA3vSC8t8nee8mxv4ccDi92d4vTLDu3w/X5/1E2t80wbEDHNXct0dU1fJR2rX3cD13TAj/OfBJevdqdZJNThQbkCVJknrOAA5OMj/JzvRC6VndsTPpfeltY0B+e/d3U04HDkhy1yR3B/4CoKp+B1ya5CD4/brnvcYqKsmuwJVV9Vl6IXfvEcfvB9xcVV8CPtgcvwG4+8Z2VfVD4IH0lmwcPdaYY3hGknt263oPAFbRuw8HJNk2yd2A53PHvXlQkn277UOBH4zS51nAn3RriOcDhwCndcfWdTPnAKcABya5d3fd9+zuTV9J5gEPrKrvA+8CdgC221R7l1hIkiT1fIPesonz6M10vrOqft0dOwP406r6WZJfAPdkjIBcVeckOabr6zf0lklsdBjw6STvARbSW8d83hh17Q+8I8k64EbgZSOOPwr4YJINwDrgdd3+lcB3kvxPVT2l2/cVYElVXTvGeGM5CzgWeADwpaoaht4XDLnjw8TnqupHSRYDFwNv6NYfXwR8emSHVXVFkmXA9+nNEp9YVd9sruH8JOdU1WHdPTu5C7zrgDcAvxhH3fOBLyXZoRvjY2M99cKALEmS5rSq2q77W8A7utfINp8HPt9trwPuNo5+3w/c6SkXVXUpveUaI/cfvom6jgKOGqX94m7zpO418vjHgY+P2L0f0PfpFc3Yl9H7Ut5Gv6qqA0Zp/8/AP4/S1e1V9ZIxat/4/mhGmdWuqnfRm/Hd+P4Y7li7PFp/V7X1VtWHmmb7jVLfqFxiIUmSNMsl2THJT4FbquqUQdcz3TmDLEmStJmS3Ive2tiRnrbxaQ/TQbec4OHtvonWXlVHAkdOYMzL+MPZ5xnDgCxJkrSZuiC5ZNB1bI6ZXPtkc4mFJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUWDDoAjS5drttw6BL6Gvhoulf401vetWgS+hr/e+m/33c/qgvDLqEvq475BWDLqGvhfeeP+gS+rr1ihp0CX3NWzj9a9x2r+0HXcKY6pbbBl2CNCmcQZYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGpWKQfEAACAASURBVAZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkjTnJTkuydlJLkyytNv3qiQ/TXJWks8m+US3f+ckxyZZ3b2eONjqtbUtGHQBA7Z2aGho+SaOLR4aGloGLJrCera69/DwQZcgSdJAdYF3abNrZVWtHNHslVV1TZK7AquTnAj8DbA3cAPwPeC8ru1HgQ9X1Q+SPAg4CdhjUi9CU2pOB+Th4eEVmzrWBedFw8PDy6esoElw3H0Pfd+ga5AkaZC6MDwyEI/05iTP77YfCLwUOK2qrgFI8lX4/azT04FHJtl47vZJtquqG7du5RqUOR2QJUmSkuxPL/TuW1U3JzkV+AmbnhWeBzy+qtZOTYWaaq5BliRJc90OwLVdON4deDxwN+BPktwjyQLghU37k4E3bXyTZMmUVqtJZ0CWJElz3XeABUl+DKwAzgQuB/4BOAtYBVwGXN+1fzMwlOT8JBcBr53yijWpXGIhSZLmtKq6FfizkfuTDFfVym4G+RvAcV37q4CDp7ZKTSVnkCVJkka3PMm5wAXApXQBWbOfM8iSJEmjqKq3D7oGDYYzyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUWDLoATa4Drzlt0CX0de3Q3oMuoa8F+z1u0CX0lQc9YtAl9HXdIa8YdAl97Xj0FwZdQl+3r/neoEvoa+GZZwy6hL6ycOGgS+jrttX/NegSxnSXvRcPugRpUjiDLEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRpzkqyY5LX92mzOMmh4+hrcZILNqOGy5LsNIH2hyf5xGaMM67rmEB/+yc5YWv1N50YkCVJ0ly2IzBmQAYWA1stWA7QYiZ4HUkWTE4p05sBWZIkzWUrgN2SnJvkg93rgiRrkhzctHlS1+at3UzsGUnO6V5PGM9ASeYn+VDX//lJ3tQcflPX15oku3ft75nkuK7tmUkePUqfOyc5Nsnq7vXEbv+fdPWem+RHSe4+ynXM7653dTfGa7pz9++u73jgoiSLknyhq+1HSZ6y2Xd7hpiTnwoGYWhoaBmwaNB1SJI01yRZCixtdq2sqpXd9jJgz6pakuSFwGuBvYCdgNVJTu/avL2qntP1ty3wjKpam+RhwNHA0DhKWUpvFndJVd2e5J7Nsauqau9uucfbgb8E/hb4UVUdkOSpwBeBJSP6/Cjw4ar6QZIHAScBe3R9vKGqViXZDlg7ynUsBa6vqsck2QZYleTkrt+9u/tyaZK/AqqqHtWF95OTPHwc1ztjGZCnzqLh4eHlUz3ogrvc/31TPaYkSdNJF4ZX9m0I+wFHV9V64MokpwGPAX43ot1C4BNJlgDrgfGGxacD/1JVt3d1XdMc+3r392zgBU09L+zafi/JvZJsP0qfj0yy8f32XSBeBfxzki8DX6+qXzVtNvpT4NFJDuze7wA8DLgNOKuqLm3q+HhXx0+S/GIC1zwjGZAlSZIm5q3AlfRmmufRm53dUrd2f9czsXw2D3h8VY2sYUWSE4Fn05sZfuYo5wZ4U1Wd9Ac7k/2BmyZQw6zjGmRJkjSX3QDcvds+Azi4W5u7M/Bk4KwRbaA303pFVW0AXgrMH+dY3wVes/GLbyOWWIzmDOCwru3+9JZhjJzNPhn4/VrmblabJLtV1Zqq+gCwGth9lOs4CXhdkoXdOQ9Pcrc+dTwceBBwcd+rncGcQZYkSXNWVV2dZFX3eLZvA+cD5wEFvLOqfp3kamB9kvOAI4FPAccmeRnwHcY/2/o5eksTzk+yDvgsMNbj2pYDRyQ5H7gZePkobd4MfLJrswA4nd466rd0X6bbAFzYXduGEdfxUXpros9Jb/3Fb4EDRhnjU8Cnk6wBbgcOr6pbR1myMWsYkCVJ0pxWVSMfffaOEcfXAU8d0aZ9osS7unaXAXuOMc7twNu6V7t/cbM9DOzfbV/DKIG1qo6kF3CpqquAg0dp86aR+zojr+N/d6/Wqd1rY19rgVeMMsYftJtNXGIhSZIkNZxBliRJ2oq6L8R9YMTuS6vq+YOoRxNnQJYkSdqKuqdCnNS3oaYtl1hIkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEkNA7IkSZLUMCBLkiRJDQOyJEmS1DAgS5IkSQ0DsiRJktQwIEuSJEmNBYMuYBpbC+zer9HQ0NAyYNE4+lu8pQVtjufvMjSIYSfkulU3DrqEvnYaumXQJfS17sgjBl1CXwvvPX/QJfR1+5rvDbqEvhY86qmDLqGvdecPD7qEvvLghw66hL4W3jzN/9uz/XaDrkCaFAbkTRgeHl4xNDS0fBxNFw0PD/dtN86+JEmSNGAusZAkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJc1aSHZO8fiv2d2qSoW77W0l23Ap9DiX52CaOXZZkpy0do+vryCQHbo2+JrPPqWBAliRJc9mOwJ0CcpIFW9pxVT27qq7bkj6SLKiq4ap685bWo/EzIEuSpLlsBbBbknOTrE5yRpLjgYuS/F2St2xsmOT9Sf5XknlJPpXkJ0m+280U32mWtJ3dTfKyJOcnOS/Jv3b7dk5ybDfu6iRP7PYvT/KvSVYB/5pk/yQndMfuleTkJBcm+RyQZryXJDmru5bPJJnfvY5MckGSNUneOp6bkmSfJKclOTvJSUl2SbJ7krOaNouTrNlU+4n/U0wfW/zpSOO2dmhoaPlUD/pg7j/VQ0qSNK0kWQosbXatrKqV3fYyYM+qWpJkf+DE7v2lSRYDXwc+kmQe8GLgscALgMXAI4F7Az8Gjhhj/D8C3gM8oaquSnLP7tBHgQ9X1Q+SPAg4CdijO/ZIYL+quqWra6P3AT+oqr9L8ufAq7ox9gAOBp5YVeuSfAo4DLgQuH9V7dm167vkI8lC4OPA86rqt0kOBt5fVa9McpckD66qS7vxjtlUe+CV/caargzIU2R4eHjFIMY9aNfnvW8Q40qSNF10YXhl34Y9Z3Xhj6q6LMnVSf4YuA/wo6q6Osl+wFeragPw6yTf79PnU7v2V3X9XtPtfzrwyOT3k8DbJ9mu2z6+qm4Zpa8n0wvoVNWJSa7t9j8N2AdY3fV3V+A3wL8DD0nycXrh/+Rx3INHAHsC3+36mg9c0R37Cr1gvKL7e3Cf9jOSAVmSJOkON414/zngcOC+jDFLvJnmAY+vqrXtzi5kjqyjnwBHVdVf3+lAshfwTOC1wIvoP7Mb4MKq2neUY8cAX03ydaCq6pIkjxqj/YzkGmRJkjSX3QDcfYzj3wCeBTyG3hIIgFXAC7u1yPcB9u8zxveAg5LcC6BZYnEy8KaNjZIsGUe9pwOHdu3/DLhHt/8U4MAk9944RpJduzXQ86rqWHrLPPYexxgXAzsn2bfra2G3TISq+i9gPfA39MLymO1nKmeQJUnSnNUtmViV5ALgFuDKEcdv65ZQXFdV67vdx9Jb0nAR8EvgHOD6Mca4MMn7gdOSrAd+RG9W+s3AJ5OcTy+TnU5vlncsfwscneRC4D+A/+7GuCjJe4CTu/XS64A3dNf0hW4fwJ1mmEep97buS4cfS7JDV9tH6K1nhl4w/iDw4HG2n3EMyJIkaU6rqkM3dawLlo8HDmrab0jy9qq6sZsVPgtY0x3bv2m3uNk+CjhqxLhX0VvDO7Ke5SPenwqc2m1fDfzpJq7jGO6Y1W2NZ9aYqjq82T6X3nrn0dp9CPjQiH2jtm/7nEkMyJIkSaNI8kjgBOAbVXXJiMMndE+EuAvwf6rq11NeoCaNAVmSJGkUVXUR8JBNHNt/aqvZepJ8EnjiiN0fraovDKKe6ciALEmSNIdU1RsGXcN051MsJEmSpIYBWZIkSWoYkCVJkqSGAVmSJElqGJAlSZKkhgFZkiRJahiQJUmSpIYBWZIkSWoYkCVJkqSGAVmSJElqGJAlSZKkxoJBF6DJtYj5gy6hr5uu32bQJfS13TdXDbqEvrIggy6hr1uvqEGX0NfCM88YdAl9rTt/eNAl9LXwsHcOuoS+bl99wqBLmPHqyqsGXYI0KZxBliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZI0pyVZnOSCzTx3/yQnbO2aNFgGZEmSJKlhQJYkSYIFSb6c5MdJvpZk2yRPS/KjJGuSHJFkG4Akz0rykyTnAC/o9s1LckmSnZv3P9v4fqQkRyb5dJIzk/y8m4k+ohv/yKbdp5MMJ7kwyd82+1ckuSjJ+Uk+1O07KMkFSc5Lcvrk3arZb8GgC5jm1g4NDS3v02bxlgwwNDS0DFi0JX2MZQ8eNFldS5I0IyRZCixtdq2sqpUjmj0CeFVVrUpyBPA24DXA06rqp0m+CLwuyb8AnwWeCvwMOAagqjYk+RJwGPAR4OnAeVX12zFKuwewL/Bc4HjgicBfAquTLKmqc4F3V9U1SeYDpyR5NHA58Hxg96qqJDt2/b0XeGZVXd7s02YwII9heHh4Rb824wjQ/SwaHh7e0j426aW7vuB9k9W3JEkzQReGRwbikX5ZVau67S8BfwNcWlU/7fYdBbwBOLXbfwlAF4o3hu8jgG/SC8ivBL7QZ8x/7wLuGuDKqlrT9XkhvQm4c4EXdQF/AbAL8EjgImAt8Plu/fPGNdCrgCOTfAX4ep+xNQaXWEiSJEGNeH/dhDuo+iVwZZKnAo8Fvt3nlFu7vxua7Y3vFyR5MPB2erPYjwZOBBZV1e1d/18DngN8pxv/tcB7gAcCZye510SvQT0GZEmSJHhQkn277UOBYWBxkod2+14KnAb8pNu/W7f/kBH9fI7eDPRXq2r9Fta0PXATcH2S+wB/BpBkO2CHqvoW8FZgr27/blX1w6p6L/BbekFZm8ElFpIkSXAx8IZu/fFFwJuBM4GvJlkArAb+papu7ZY8nJjkZuAM4O5NP8fTW1rRb3lFX1V1XpIf0Qvlv6S3hIJuvG8mWQSE3nppgA8meVi37xTgvC2tYa4yIEuSpDmtqi4Ddh/l0CnAH4/S/jubaA+92dzzquonfcY8fMT4e27i2OGM7rGj9PmCscbU+BmQJUmStoIky4DX0XuShWYwA7IkSdJWUFUrgD94AlaSdwMHjWj61ap6/5QVpgkzIEuSJE2SLggbhmcYn2IhSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1Fgy6gNlmaGhoGbBoAqcsnqRSALiNDZPZ/VZxy9qFgy6hr/U31qBL6GvdDRl0CX3NWzj972MWTv//PebBDx10CX3dvvqEQZfQ14LHPGfQJfS17uILB13CmOqa6wddgjQpDMhb36Lh4eHl4208NDQ07raSJEmafC6xkCRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJEmSGgZkSZIkqWFAliRJkhoGZEmSJKlhQJYkSZIaBmRJkiSpYUCWJElzVpIdk7x+CsY5IMkjJ3scbR0GZEmSNJftCIw7IKdnc/LTAYABeYYwIEuSpLlsBbBbknOTfDjJKUnOSbImyfMAkixOcnGSLwIXAA9M8jfdvh8kOTrJ27u2uyX5TpKzk5yRZPckTwCeC3ywG2e3gV2txmXBoAsQa4eGhpZPVue78YDJ6lqSpBkhyVJgabNrZVWt7LaXAXtW1ZIkC4Btq+p3SXYCzkxyfNfuYcDLq+rMJI8BXgjsBSwEzgHO3tg38NqquiTJ44BPVdVTu35OqKqvTerFaqswIA/Y8PDwisns/+BdD3jfZPYvSdJ014XhlX0bQoB/SPJkYANwf+A+3bFfVNWZ3fYTgW9W1VpgbZJ/B0iyHfAE4KtJNva5zda5Ck0lA7IkSVLPYcDOwD5VtS7JZcCi7thN4zh/HnBdVS2ZpPo0RVyDLEmS5rIbgLt32zsAv+nC8VOAXTdxzirgL5Is6maNnwNQVb8DLk1yEPz+C317jTKOpjkDsiRJmrOq6mpgVZILgCXAUJI1wMuAn2zinNXA8cD5wLeBNcD13eHDgFclOQ+4EHhet//fgHck+ZFf0pv+XGIhSZLmtKo6dBzN9hzx/kNVtTzJtsDpdF/Sq6pLgWeNMsYqfMzbjGFAliRJmriV3Q9/LAKOqqpzBl2Qth4DsiRJ0gSNc9ZZM5RrkCVJkqSGAVmSJElqGJAlSZKkhgFZkiRJahiQJUmSpIYBWZIkSWoYkCVJkqSGAVmSJElqGJAlSZKkhgFZkiRJahiQJUmSpIYBWZIkSWosGHQBs8DaoaGh5c37xQOqY1T7bdhu0CX0tevQLwddQl+Zl0GX0Nf6W2rQJfS17V7bD7qEvm5b/V+DLqGvhTffMugSZoV1F1846BL6WviSdw26hDHd9ol3D7oEaVIYkLfQ8PDwivb9iLAsSZKkGcYlFpIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRJUsOALEmSJDUMyJIkSVLDgCxJkiQ1DMiSJElSw4AsSZIkNQzIkiRpTkuyOMkFo+w/NcnQVuj/8CSf2NJ+NHUMyJIkSdNckgWzcazpyoAsSZIEC5J8OcmPk3wtybbtwSSHJP+/vTuPs6Oq0z/+eRJAFAiguKHsCogIyCKgjIr7IKADIgNug8wIyijuwKCCCyruiOOCG4vID1BRwAUE2QxCSAgQFtER3FFEIUT25fn9ceom1Z1ebiDpU7f7eb9e/UpX1e3uJ919+37r1Knv0TxJV0s6so/9+0j6laRZwHPG+sKSjpX0ZUmzm4/Zudn/H5JOl/Qz4FxJK0n6hqRZkuZKekXzuKc3+66QdJWkpzaP/aGkK5tsezaP/a2kNZr3t5Z0fvP+4ZJOkDQTOEHSYyV9V9JlzduY/4fJZsqfISwLW2+99cHAirVzALyBjWtHiIiIqErSm4A3tXYdY/uYYQ/bCNjX9kxJ3wDe0vr4NYEjga2AW4GzJb0SmDXK/kuBDzb75wPnAXPHibku8CxgA+A8SU9p9m8JbGb7H5I+CvzM9hslrQbMknQOsD9wlO0TJa0ATAd2Av5s++XN/2HVPr5VmwA72L5L0reBz9r+uaS1gbOAp/XxOSaFFMjLxoqzZ88+vHYIgKPXeu1htTNERETU1BTDwwvi4f5ge2bz/reAt7WObQOcb/tvAJJOBJ4LeJT9DNt/MrDhOF//FNsPAr+WdAMsHOH6qe1/NO+/BNhV0rub7RWBtYFfAIdKejLwPdu/ljQP+HQzqn2m7YvG+foAp9u+q3n/RcAmknrHZkha2fY/+/g8Ay8FckREREQpdsfarvX172jtE7C77euHPfY6SZcCLwd+JGk/2z+TtCVlJPkjks61/SHgfhZNsR1+tbv9taYB29m++yH+fwZa5iBHREREwNqStm/e3xv4eevYLOB5ktaQNB3YC7hgjP2XNvsfI2l5YI8+vv4ekqZJ2gBYHxheBEOZ5vBWNcO6kp7Z/Ls+cIPtzwM/ADZrpoXcaftbwCcpUzUAfkuZ+gGw+xh5zgbe2tuQtEUf/4dJIwVyRERERClID5B0HbA68KXeAds3AQdT5hJfCcyx/YNx9h9OmfowE7iuj6//e0rB/WNg/1FGbj8MLA9cJemaZhvg1cDVkq4ANgWOB55BmaN8BXAY8JHmsR8EjpI0G3hgjDxvA7Zubvq7ljLPecrIFIuIiIiY0mz/Fka8q/35rcecBJw0wseOtv+bwDeXIMY5tocUobaPBY5tbd8F7DfC1/o48PFhu89q3oY/9iJGmA9t+/Bh27cAe/YbfrLJCHJEREREREtGkCMiIiImgKRDWXw+8qm2/6NCnBhDCuSIiIiICWD7COCI2jlifJliERERERHRkgI5IiIiIqIlBXJEREREREsK5IiIiIiIlhTIEREREREtKZAjIiIiIlpSIEdEREREtKRAjoiIiIhoSYEcEREREdGSAjkiIiIioiUFckREREREy3K1A0xCdwMb1w7R88fpD9SOMK4H76ydYHyP2Gjl2hHGNX21e2pHGJfvurd2hHGtsOW6tSOMb0b3fx/911tqRxiX/zG/doRx3fuFQ2tHGNMK/31E7QgRy0RGkJey2bNnfxz4be0cEREREfHQpECOiIiIiGhJgRwRERER0ZICOSIiIiKiJQVyRERERERLCuSIiIiIiJYUyBERERERLSmQIyIiIiJaUiBHRERERLSkQI6IiIiIaEmBHBERERHRkgI5IiIiIqIlBXJEREREREsK5IiIiIiIlhTIEREREREtKZAjIiIiIlpSIEdEREREtKRAjoiIiClN0rqSrh5h//mStl4Kn/8/JH3h4X6eZUnS8yWd2fXPOVFSIEdERERMEZKWq51hEKRAjoiIiIDlJJ0o6TpJ35H0qPZBSXtJmifpaklH9rF/H0m/kjQLeM5YX1jSHs3HXynpwmbfkFFnSWdKen7z/r69zy3pq73HSdpF0qWS5ko6R9Ljm/2HSzpB0kzghPG+EZJWkvSN5vPPlfSKZv8lkp7eetz5krYe7fGDLGcRy8bdW2+99eG1QwC8kKfWjhAREVGVpDcBb2rtOsb2McMethGwr+2Zkr4BvKX18WsCRwJbAbcCZ0t6JTBrlP2XAh9s9s8HzgPmjhHxA8BLbf9J0mrj/F/WBN4PbAksAH4GXNkc/jmwnW1L+k/gvcC7mmObADvYvmusz984FPiZ7Tc2eWZJOgc4GXg1cJikJwJPtD1b0kdHefzASoG8DMyePfvjtTP0HLTuXofVzhAREVFTUwwPL4iH+4Ptmc373wLe1jq2DXC+7b8BSDoReC7gUfYzbP/JwIZjfO2ZwLGSTgG+N07OZwEX2P5H87lPbX3uJwMnN8XrCsCNrY87vc/iGOAlwK6S3t1srwisDZwCnA0cRimUvzPO4wdWCuSIiIiIUuyOtb3svrC9v6RtgZcDcyRtBdzP0KmwK/bxqY4GPmP79GY6xuGtY3csQSQBu9u+frED0t8lbQbsCew/1uN7UzwGUeYgR0RERMDakrZv3t+bMl2hZxbwPElrSJoO7AVcMMb+S5v9j5G0PLDHWF9Y0ga2L7X9AeBvwFrAb4EtJE2TtBZl5BjgsuZzr97ccLd761OtCvypef8ND+F70HMW8FZJavI9s3XsZMrUjVVtX9XH4wdSCuSIiIgIuB44QNJ1wOrAl3oHbN8EHEyZS3wlMMf2D8bZfzjwC8r0ievG+dqf7N3oB1zcfK6ZlCkS1wKfBy5vsvwJ+CilOJ9JKaTnN5/ncOBUSXOAWx7qNwL4MLA8cJWka5rtnu8A/06ZbtHP4weS7Am7ghAVHLTuXp3/Ab9noz/XjjCuR2y0Su0I4/Id99SOML7lun9Ovtz6a9aOML4ZK9dOMC7/9eG8Nk+QBzv/5xGmqXaCMa3w30fUjtCX5ddYv9vfyCUkaWXb/2xGkE8DvmH7tNq5JpPuv1pFRERERNvhkq4ArqaMMn+/cp5JJzfpRUREREwASYey+HzkU20v0VC87XeP/6hRM7yU0pqu7Ubb//ZQP+dklAI5IiIiYgI0hXDVeSm2z6LcVBdjyBSLiIiIiIiWFMgRERERES0pkCMiIiIiWlIgR0RERES0pECOiIiIiGhJgRwRERER0ZICOSIiIiKiJQVyRERERERLCuSIiIiIiJYUyBERERERLSmQIyIiIiJaZLt2hhgwkt5k+5jaOcaSjEtHMi4dybh0JOPD1/V8MBgZY/LLCHI8FG+qHaAPybh0JOPSkYxLRzI+fF3PB4ORMSa5FMgRERERES0pkCMiIiIiWlIgx0MxCHPDknHpSMalIxmXjmR8+LqeDwYjY0xyuUkvIiIiIqIlI8gRERERES0pkCMiIiIiWlIgR0RERES0pECOiIh4iCSt18++iBgsuUkv+iLpQOCbwALga8AzgYNtn101WEPSvra/3tqeDrzP9gcrxupl2XKs47Yvn6gsS0rSNGBl27fXzjKcpEfZvrN2jkEjabexjtv+3kRlWVKSVgfWsn1V7Sw9ki63veWwfXNsb1Ur01jyvInoz3K1A8TAeKPtoyS9FFgdeB1wAtCJAhl4oaTdgX2BRwPHAhdUTbTIp8c4ZuAFExWkH5K+DewPPABcBsyQdJTtT9ZNVkh6NuUkbWVgbUmbA/vZfkvdZCDpDMrPdES2d53AOKPZZYxjBjpVIEs6H9iV8no1B7hZ0kzb76yca2Pg6cCqw046ZgAr1kk1ujxvIpZMCuTol5p/dwJOsH2NJI31ARPJ9t6S9gTmAXcAe9ueWTkWALZ3rJ1hCW1i+3ZJrwF+DBxMKUw6USADnwVeCpwOYPtKSc+tG2mhT9UOMB7b+9TOsIRWbX4f/xM43vZhkrowgrwRsDOwGkNPOhYA/1Ul0djyvIlYAimQo19zJJ0NrAccImkV4MHKmRaS9FTgQOC7wNOA10ma26VLiZKWB94M9F6Uzge+Yvu+aqFGtnyT9ZXAF2zfJ6lTc7Fs/2HY+dkDtbK02V541ULSCsCGzeb1U+heNgAAIABJREFUXfs5S3o88FFgTdv/KmkTYPv2VKWOWE7SE4FXA4fWDtNj+wfADyRtb/sXtfP0Y0CeN48E1rZ9fcVIEblJL/q2L2UkcZum6FwB6NJI1BnAB2zvBzwP+DVlekCXfAnYCvhi87ZVs69rvgL8FlgJuFDSOkCX5iD/oblcbEnLS3o3cF3tUG2Snk/5Hfxfys/6Vx0ares5FjgLWLPZ/hXw9mppRvchSs7/s32ZpPUp39uu+LukcyVdDSBpM0nvqx1qBIPwvNkFuAL4SbO9haTT66aKqSo36UXfJD0JWIfWlQfbF9ZLtIikGcNvJJO0oe1f1co0nKQrbW8+3r4ukrSc7ftr5wCQtAZwFPAiytSfs4EDbf+9arAWSXMo03yub7Y3BE7q0o1bki6zvU1zpeWZzb4rbG9RO9sgkXQB8B7K1aDe9/Fq25vWTTbUAD1vXgCc3/pezrP9jLrJYirKFIvoi6QjgT2Ba1l0Wc5AJwpk4JGSPgs8yfbLepeLKaNiXfGApA1s/wagGQnrxCXONkmPAHYH1mXo34gPVQk0jO1bgNfUzjGO5duXiG3/qpm20iV3SHoMzc1RkrYD5teNtDhJnwA+AtxFGVncDHiH7W9VDbbIo2zPGjZ1oRMnkz1NV5+jbHf9eXOf7fnDvpcZxYsqUiBHv14JbGT7ntpBRnEspQ1db47ir4CTgS7Np3w3cJ6kGygjOOvQrWkqPT+gFEpzgM79vCV9foTd84HZzbzQLpgj6WtAr4h7DTC7Yp6RvItyw9YGkmYCjwVeVTfSiF5i+72S/o0y9Wc3yol5VwrkWyRtwKITjVcBN9WNNJTtByStI2kF2/fWzjOGayTtDUxv7it5G3Bx5UwxRaVAjn7dACxPBwumxhq2T5F0CIDt+yV1ZnS2GcHZHHgq5e53KDdudfH7+WTbL6sdYgwrAhsDpzbbuwM3AptL2tF2F+bR7g8cQHmBB7iIMhe5M2zPkfQ8yu+j6OCNhI3e69TLgVNHGGGs7QDgGGBjSX+i/C6+tm6kEd0AzGzm9N7R22n7M/UiLeatlEGOe4BvU+aef6RqopiyUiBHv+4ErpB0Lq0i2fbbRv+QCdXpy8XNCM5etj8LdKFF1VgulvQM2/NqBxnFZsBzbD8AIOlLlAJ0B0qbv6qak6ErbW8MdKn4GKJplfb/gJN703466kxJv6RMsXizpMcCd1fOtJDtG4AXSVoJmGZ7Qe1Mo/hN8zYNWKVyltFsbPtQOtStJKau3KQXfZH0hpH22z5uorOMpFmt7mhgU+BqmsvFHVtx67OUUfiTGTqC06mV9CRdCzyFMhJ2D2V00bY3qxqsIel64Fm25zfbqwKzbG/UvuGsJkk/AN5q+/e1s4ym6U6yZ/P2IOX38pQuZpb0aGB+c6L5KGCG7b/UzgUgaaQFS+YDc2xfMdF5Bpmk84AnAN+hnLhdXTlSTGEpkGNSkLQH5XLcWpRL7tsC7+9S8dn88R/Otru2kt46I+23/buJzjISSfsC76P0kRalr/RHgZOAw22/p166QtKFlOXYZzH0ZKiTK4I18z3fD7zG9vTaeYaTtCmwCa0V6mwfXy/RIiorT25NaTUJZfGQqyg3uZ5q+xOVog3R/P1Z7AW/g39/nkDpeb0nZVXCk21nmkVMuBTIMSZJp9h+taR5jPzHtSujilfZ3kzSDsCHKSszfcD2tpWjDYxeq7xmtG4xtv8x0ZlGI2lNynLn11GWzv1jV1oOAjRzexfTXhChC4aNIj9AKUbGWhp9wkk6DHg+pUD+EfCvwM9td+KGwuZkaCfb/2y2VwZ+CLyMMoq8Sc18PZLaLQZXpAwk3G/7vZUijUnSM4D3AnvaXqF2nph6Mgc5xnNg8+/OVVOMr3dD3suBr9r+oaROjToMwMpl36b8nOdQTobad0IZWL9GqOFUlhw+EHgyZVGB7YBfUPqndoLtC5ri86m2z2mmBXRqZFbSpZQpP6cCezRzabvoVZQbXOfa3qd5HnWlgwXA4xh68/J9wONt3yWpMzfh2p4zbNdMSbOqhBmFpKdRTtZeBdxCmfbzrqqhYspKgRxjsn1T8+/vhr3gP5Ju/f78SdJXgBcDRza9fLu2UuSxdLgVne2dm3/Xq51lHAcC2wCX2N5R0saUE4/OkPRfwJuARwMbAE8Cvgy8sGauYV4/IMv53mX7QUn3S5oB3EyZStUVJwKXNvPOAXYBvt3ctHdtvVhDDbsyNI2ykueqleKM5huUG0dfYvvPtcPE1NalAic6bIQX/CfTrRf8V1MuaX7K9m2SnkhZ3apLOt2Krkelh9ZrgPVsf1jS2sATbHdltOlu23dLQtIjbP9S0kbjf9iEOgB4FnApgO1fS3pc3UiLuU3S1+nuFY2e2ZJWA75KubrxT8oVg+qa58qxwI+B5zS797fd63ndpYU52leG7qfchLtv1UTD2N6+GXxZu3aWiBTI0a9Ov+DbvhP4Xmv7JjrWrJ+Ot6Jr+SKlq8ELKPO5FwDfpYzadsEfm4Lp+8BPJd0KdOIGwpZ7bN/b69craTm6tyLYsXT4ikaP7bc0735Z0k8oHSw60Z3GtiX9qFkKuWsLwQwxAFeGkLQL5f6RFYD1JG0BfKirN7fG5Na1S9DRXfe0V2Dq6At+172ToSuXHc+ihSS6ZFvbB9D0mrV9K+UFqxNs/5vt22wfTum88HXKSo9dcoGk/6Esgf5iyjzfM8b5mIm2hu1TKCdD2L6fDi59DiBpM0m7AlsCT5G0W+1MLZdL6srJ46gk7SFpleb990n6XtMes0sOpwzE3AbQtMnrfGEfk1NGkKNfw1/w30L3XvC77hpgyMpldPMk9b5msYveSPdjaYqorulaV4iWgymXr+cB+wE/sv3VupEWMxBXNCR9g7I4zDUs+j00rStGlW0LvEbS7ygt/TrVN7zl/bZPbTr9vAj4JPAlSv6uuG+ElRIzEBNVpECOfg3CC37X/cL2lpQXegAkXU4ZFeuSzwOnAY+TdATljvL31400cN5q+yjKvFkAJB3Y7OuK4Vc0Hkv5WXfNdl1plTaKl9YO0Kd2p59jutjpB7hG0t7A9KY399uAiytniikqfZCjLyO9uHfwBb+Tmsb3T6K0ptqbRe3TZgBfbpYk7pSmM8QLKVnPtX1d5UgDRdLlzclQe18nVvlra6ZKLbyiYfu+ypEW09xI+GnbnekIMZLmnoz2QiadWpFQ0pnAnyidfrakLN09y/bmVYO1NO0QDwVeQvmdPAv4sO3OLC0eU0cK5OjLoLzgd5HKMt3/QVlt6zIWFcgLgGNtd+VSMQCSTrD9uvH2xeIk7UU5CdoBuKh1aAbwgO3qXV/Gm7/bwd/H51FGuv9CN5c+3xX4NLAmpQXdOsB1tp9eNdgwTfH5MmBec5P1E4Fn2D67crSITsoUixhT6wV/PUmntw7NADqzslqX2T4OOE7S7ra/WztPH4a8sDfzkbca5bEx1MWU7ilrUIqmngWU5Ye7YJcxjnVpbm/P1ymrJs6jm3PhP0xZrOYc28+UtCPw2sqZRvJE4Ie275H0fMq87q4s130GY8w1TheLqCEFcoxnEF7wB8WTm4UOFlDmpm4JHNyVEZymP3PvRszbWTTSfS9wTLVgA8T274DfSXoRixa42BDYmFLgVWd7n34eJ+kNzcldbX+zffr4D6vmPtt/lzRN0jTb50n6XO1QI/gusLWkp1Cezz+grJ65U9VUxadqB4gYLlMsom/NEq+9dkazbN9cM8+gkXSl7c0lvRTYH3gfcMLwqSu1SfqY7UNq5xhkkuYA/wKsDsykTK2513aXFo4Y00jTqirl+CKwGqVrzsKlm7syFUTSOZQ2gx+jDCTcDGxt+zljfuAE6/08Jb2XcvJ29KBNk5P0Xdu7184RU0NGkKMvkvagnOWfTxlZPFrSe2x/p2qwwdIbkd0JON72NRrWz6gLbB/SzKt8brPrfNtn1sw0gGT7Tkn7Al+0/QlJV9QOtYS68rv5SEph/JLWvi5NBbkSuBN4B2XlvFWBlasmGtl9zZS517Noms3yFfM8FOvXDhBTRwrk6Nf7gG16o8ZNb9xzgBTI/Zsj6WxK4/tDmqb9nZtTKeljlGb9Jza7DpT0bNv/UzHWoJGk7SkFU2853+kV8zwUnbi8ON6UEEmH2P7YROUZwY62H6Q8l49rMnVx+tk+lCtXR9i+UdJ6wAmVMy2pTvxOxtSQKRbRF0nzmuVUe9vTgCvb+2JszfdsC+AG27c1izQ8qSvL5vY0L+5bNC/6vZv05nala8AgaDovvAuYaftISesDb7fdxZUTRzQol99rTQWR9GbKgkkbAP/XOrQK5efeuRv1JD0SWNv29bWzPBRdmfYTU0NGkKNfP5F0FnBSs70n8KOKeQaRgU2AnYEPASvR6pvaMauxqEvJqjWDDKJmhb8LmtZa2L6Bbi4rPpaZtQP0qdZUkG8DP6bMPT64tX+B7c51+JG0C2Wa3AqUrkRbAB8asA4RXZn2E1NAF5e5jY5p5sl+HvgKpTXQZpSVmA6qGmzwfBHYHtir2V4A/G+9OKP6KDBX0rGSjgPmAEdUzjRQJG0v6Vrgl8325s3NZp0h6UBJM1R8XdLlkhbO87X93zXzLYEql0Ftz7f9W9t72f5d661zxXHjcMrUqdsAbF9Bh+b0Spou6cRxHpbXnJgwGUGOcdm2pB810ym6cmPMINq2uYt8LoDtWyWtUDtUWzMN5EFKX9dex5KDbP+lXqqB9DnKEsSnA9i+UtJzx/6QCfdG20c1XVVWp/QaPgHoRNvBJZBRxf7cZ3v+sPuCO3MPhO0HJK0jaQXb947ymEH73YwBlgI5+nW5pG1sX1Y7yAC7r5nPa1h4o2NnXqAAmr6977V9Ck1xFw+N7T8MK0YeqJVlFO2uKid0tauKpMfY/vsYDzl1wsIMtmsk7Q1Ml/RUypSfiytnGu4GYGazKNUdvZ22P1MvUkxVmWIR/doWuETSbyRdJWleR+/U7rLPA6cBj5N0BPBzynSGrjlH0rslrSXp0b232qEGzB8kPRuwpOUlvRu4rnaoYXpdVXYCzupqVxXK351TJe00UgFvu4vPoS56K2WVzHso86fnA2+vmmhxvwHOpNQmq7TeIiZculhEXyStM9L+ZuWw6JOkjYEXUkbvzrXdtaIJSTeOsNu2OzNfseskrQEcBbyI8rM+GzhwnJHQCTVAXVVE+T6+kTLt5xTgWNu/qhpsgDRXrs6xvWPtLBGDIgVy9E3SlsAOlCkCM21fXjnSQBhv9LXDN/XEJNYUnq8B1rf9IUlrA0+wPatytFFJ2hH4FqUDzJWUpdp/UTfVYJB0LrCb7fm1s4ymmXb2XspI98IOP7ZfUC1UTFmZgxx9kfQBYA8W3aT3TUmn2v5IxViDYg7lpGKk+Z2mQ3eSA0hakdLftXcydBHwZdt3Vw02ACQdzRhdFTrWB/mLlCkVL6C0HVwAfJdFN2d2QjOy/VrKTYR/pUwVOJ0y+n0qZeGdGN8/gXmSfsrQ+b1d+p08ETiZ0gpzf+ANwN+qJoopKyPI0RdJ1wOb94qkpuH8FbY3qpssljZJp1CKpW81u/YGVrO9R71Ug0HSG8Y6bvu4icoynt6iC+0FQSRdaXvz2tnaJP2K0l3jm7b/OOzYQbaPrJNssIz2u9mx38k5treSdFVvYSJJl9nu1ElbTA0ZQY5+/Zlyyas3ivgI4E/14gwmSbsCvXZf59s+s2aeUWxqe5PW9nlNT98Yx/BiQ9LKzf5/1kk0ps53VWls5FFGclIc98/2cU1byY0pP/PrR2unVtF9zb83SXo55XUnNwhHFSmQo1/zKW2Cfkr54/piYJakz0PnLtN1kqSPUy5f95rhHyjp2bb/p2KskVwuaTvblwBI2haYXTnTQJG0KWXU89FlU38DXm/7mrrJhhjeVeVVwPvqRlpE0hksKt4XOz5gK8BVJ2knymJPv6FM91pP0n62f1w32RAfkbQqZZn2o4EZwDvqRoqpKlMsoi+DdOm4q5q2eFvYfrDZng7M7V1K7ApJ1wEbAb9vdq0NXA/cT+lm0am8XSTpYuBQ2+c1288HPmr72VWDDdPlriqSnte8uxvwBBZN+dkL+KvtFE5LQNIvgZ1t/1+zvQHwQ9sb100W0U0ZQY6+tC7Pbdjsut72fWN9TIxoNaDXtWLVmkHG8LLaASaBlXrFMYDt8yWtVDNQz7CuKjcDJ7WPdaWriu0LACR92vbWrUNnSMoVjSW3oFccN26g3GvQGb0rksPMB2bb/sFE54mpLQVy9KUZATsO+C1ltGktSW+wfWHNXAPmY8BcSedRvofPBQ6uG2lxtn8naXPgX5pdF9m+smamAXSDpPdTpllA6cJwQ8U8be2uKmsDtzbvr0a5atC1rhArSVrf9g0AktajtHmLJTNb0o8ofaRN6Up0maTdAGx/b6wPniArUuZI91ZH3B24Edhc0o62u7awSUximWIRfZE0B9jb9vXN9obASba3qptssEh6IovaaM2y/ZeaeUYi6UDgv1jU0u/fgGNsH10v1WCRtDrwQUqrPCit8g63fWu9VENJ+ipwmu0fNdv/CrzS9n51kw0l6WXAMZQTDAHrAG+yfXbVYANG0jfHOGzbb5ywMKOQdAnwHNsPNNvLUZ47OwDzht08HLFMpUCOvrTb7oy1L8Ym6UmUF/iFV2+6NgrfzJXe3vYdzfZKwC/ys55cJM2z/Yzx9nWBpEdQRhYBfmn7npp5JiNJh9j+WOUM1wPP6i1m0tywN8v2Ru12hBETIVMsol+zJX2NRTfKvIZ0Nlgiko4E9gSuYVE7LQOdKpApo3QPtLYfYORFTmIUzRWWdwPrMvRkqEsrgv1Z0vsY+pz+c8U8I5K0PLAfrfaIkr6SeyCWuj0o08Bq+gRwhaTzWTQN7aPNSfo5NYPF1JMR5OhLM4JzAEMvGX8xIzn9a0ZHNuv690zSOykrWJ1GeZF6BXCs7c9VDTZAJF0JfJky33fhyYbtOdVCDdPcrHcYiwrPC4EPduUmvZ7mxHx5yj0QUFbUe8D2f9ZLNfl0ZYS2mYb2rGbzMtudO2mLqSEFcvSt6WLxNMroZxebzHeapB8De3R00YghJG3JoqWmf257buVIA6W3IljtHJPBSKv7dXHFv0HXW1mxcobnjrS/a9PQYmrIFIvoS7Oq0ZfpdpP5TpJ0NKXQvJNy+fBcYOEococXWRGLuh1EH1ot1M6Q9BbKKHz7Z119dFbS52y/vb0QR1sHF+B4QNIGtn8DIGl9hk4BiqWjC8/z97TeX5EykjwH6NLUpJgiUiBHvz4N7Di8yTyQAnl8vbnac4DTawbph6QPUOYjfpfyovlNSafa/kjdZAOh3UINhr7gG1h/whMtrtd67lNVU/TvPZTlzntt8tYF9qkXZ9I6dfyHLFu2d2lvS1oLyNSuqCJTLKIvki6zvU1rW5S7i7cZ48OipbnR5O5WC6PpwCNs31k32VDNXOnNbd/dbD8SuML2RnWTxdLS/O4db/s1tbOMR9KKlKWHXwjcBlwGfLb3+xlja13BGlGHr2D1XmeuSXu3qCEjyNGvQWgy33XnAi8CenOQHwmcDXRq+WFKJ4MVgV4B8gjgT/XiDB5JBwAn2r6t2V4d2Mv2F+smK2w/IGkdSSsMwL0ExwO3Ax9utvemjILvUS3RYOldwXoOsAlwcrO9B3BtlUSjGFbMTwO2AC6vlyimsowgR18Gocl810m6wvYW4+2rTdL3KYuZ/JTyYvViYBbwR+j2iFNXjPKz7kSXgB5Jx1Nuuj0duKO33/ZnqoUagaRrh48gjrQvxtYswrGD7fub7eUpq2RuVzfZIpLe0Nq8H/it7Zm18sTUlhHk6IvtzPl7+O6QtKXtywEkbQXcVTnTSE5r3nrOr5RjkE2XJDcjEM2UhhUqZxruN83bNGCVylnGcrmk7WxfAiBpW9KD/aFYHZgB9G4UXbnZ1xm2e638eldd1qoYJ6a4FMgxJknvtf2J0eaxZTRxibwdOFXSnyk3cT2BsnBIp7RfpOIh+wlwsqSvNNv7Nfs6w/YHASQ9qmvz4KGs6kf5m7M8cLGk3zfb6wC/rJltQH0cmCvpPBYtwnF41UTDNAuE7EqpTeYAN0u62PY7qgaLKSlTLGJMknaxfcawS18LpZhaMs1lzd7Nbte3VwOT9GLbP62TDCSdYvvVrcJkiCw13T9J04A3UeacQ5mu8rXeDZpdIGl74OvAyrbXlrQ5sJ/tt1SOBoCkdcY6bvt3E5VlspD0BGBbyvN7lu2/VI40RG8akqT/BNayfZikq/K3J2pIgRzREbUb9Ut6ou2bRitMUpAsPZK+a3v3yhkuBV4FnN6bGy3patub1swVy46kXVm0cuIFts+omWe45uT8JZRVEw+1fVkK5KglUyxiTKMtJtDTwUUFBlnVRv22b2r+TSG87HWhHzK2/1A6aS3UmRHuWLokfZxy8+2Jza63Sdre9v9UjDXch4CzKKt3XtYsCvPrypliikqBHOMZlMUEJoOql3MkLRglgyidSmZMcKTJrAuX7v4g6dmAm6k/BwLXVc4Uy85OwBa2HwSQdBwwF+hMgWz7VFoLlti+AVh4pUXSIbY/ViNbTD0pkGNMti/o53FduGQcD4/tvjoZSFrd9q3LOk8sc/sDRwFPovS5Phs4oGqiWNZWY1EXi1VrBnmI9gBSIMeESIEcS0snLhkPuN/WDtCnc4Fqc6UniarTaQBs3wJ0fiW9WGo+xuJdLA6uG2mJVX/exNSRAjmWli5cMu60phfuy4F1aT33egsz2N6tTrIllheph++g2gGaS+wHDlvt79NZ9Gdysn1S00Ztm2bXQV3rYtGHvM7EhEmBHDFxzqAs3zwPeLBylocjL1KjGK1FHovmcW9GeefsCQ02ss16xTGA7VsldWalv1gmtmFRFwtT/iYNkpycx4RJgRxLS/5wje/JaVc06e1cO8ASmNaeTy7p0eQ1YdIakC4W4zl1/IdELB3pgxx9kbQL8MPeHdAjHH9JR0bFOkvSkcC5g/596jXzr50jHh5Jr6d0MDiVcoL7KuAI2ydUDRbLhKSrGNrFYjowtwsn7aOt1NqTFVujhmm1A8TA2BP4taRPSNp4+MFBL/omyCXAaZLuknS7pAWSbq8dajhJn5b09DEe8sIJCzOgJG0n6TJJ/5R0r6QHuvaztn08sBvwV+AvwG4pjie91Vrvd6mLxWzK0tIrUm4A/nXztgWwQsVcMYVlBDn6JmkGsBewD+Vs/5vASbYXVA02ICTdCLwCmOcOP/GaZV73oVxu7/2M59dNNVgkzQb+nTI6uzXwemBD24dUDUZ5Htu+vZlSsRjb/xhpfww2SXsBHweGdLGwfXLVYC2SLgF2sH1/s708cJHt7eomi6koBXIsEUmPAV4HvJ2yqMBTgM/bPrpqsAEg6ULg+aNNU+kaSRtRCuW9gJnAV22fVzfVYJA02/bW7WVyuzI1RdKZtnduTtjaLwC9GwnTsnGSkvREFnWxmNW1LhaSrge2752kNZ1VLrG9Ud1kMRXlhozoi6RXAP9BKYiPB55l+2ZJjwKuBVIgj+8G4HxJPwbu6e3stXnrkmZ+4sbN2y3AlcA7Je1n+9+rhhsMd0paAbhC0ieAm+jIlDbbOzf/rlc7Syx7kob3LP9j8++akta0fflEZxrDx1m8V/PhVRPFlJUR5OhL0zP167YvHOHYC22fWyHWQJF02Ej7bX9worOMRdJnKd0Yfkb5mc9qHbs+oznjk7QOZW7vCsA7KPM9/9f2b6oGY8SCaYiOFUzxMDXFZs9IVwxeMMGRxiTpCcC2lKydG+WOqSMFcvRF0pG2DxpvXww+SfsAp9i+Y4Rjq2Y+8vgkHWj7qPH21TCsYBqucwVTLB2SHgm8BdiBUnxeBHzJ9t1Vgw0jaVcW9Wq+wPag9WqOSSIFcvRF0uW2txy2b+H8yhhfU5gs9oTrWkEi6VzbLxxvX4xulOdLJ+Ygx9Qk6RTgdhb1Qd4bWNX2q+ulGmqEXs17AZcNWK/mmCQyBznGJOnNlFGHDZo+mj2rUG7civ69u/X+isDuwP2VsixG0orAo4A1mptjeou/zACeVC3YAGk6BewNrCfp9NahGUCnukNIGmlp8/mULis3T3SeWOY2tb1Ja/s8SddWSzOynRjaq/k4YC6lX3fEhEqBHOP5NvBj4GPAwa39C9IOasnYnjNs10xJs0Z8cB37UbqTrEnpSdorkG8HvlAr1IC5mHJD3hrAp1v7FwBXjfgR9ewLbE9p+wXwfMrPfT1JH0pP5Enncknb2b4EQNK2lP7DXbMai04mu9SrOaaYTLGIMaVn6tIz7Hs4DdiK0iKvUze9SXpr2vY9fJIez9CWWp0alZV0FvB6239tth9P6VCzF3Ch7U1r5oulQ9I8ytSu5YGNgN832+sAvxw2qlzVIPRqjqkjI8gxnm9TOhrMofxRVeuYgfRM7V/7e3g/cCNlFK9r/iJpFdsLJL2PsrLVR9LdoH+S9gA+BZxP+XkfLek9tr9TNdhQa/WK48bNzb5/SLqvVqhY6nauHaBftk+SdD6LTiwPSheLqCUjyBETQNI0SgP8zs/b7t18KWkH4CPAJ4EP2N62crSBIelK4MW9UWNJjwXOsb153WSLSPoisDZltT8oc+L/CLwHONP2jrWyxdSVLhbRFZ1oXB/dJ+k5klZq3n+tpM9IWrt2rkHR3HQyKPN4H2j+fTlwjO0fUvr5Rv+mDZtS8Xe69/f2AMpS4ls0b8cDB9i+I8Vx1NB0sTiQsvjUtcDbJH20bqqYqjKCHH1pOlhsDmwGHAt8DXi17efVzDVIJH0K+AXwPXf4iSfpTOBPwIsp0yvuosyh7czoZ9c1q+dtDpzU7NoTuKpLfcMlvRX4lu1ba2eJgIWvM+0uFtOBuWknGjV0bUQjuuv+pqh7BfAF2/9LafUW/duPcjn7Xkm3S1og6fbaoUbwauAs4KW2bwMeTbnsHv0z8BXKCeX2XDipAAAIuUlEQVRmwDF144zo8cBlkk6R9DJJGvcjIpa91Vrvp4tFVJMR5OiLpAuAnwD7UOaH3QxcafsZVYPFUpOOJUvPoCys0xTFL6E8r7cGTqEsL159SeyYetLFIrokXSyiX3tSFkDY1/ZfmvnHn6ycaaA0xchrgPVsf1jSWsATbXelF3I6ljxMrYV11h+EhXVsW9JfgL9QOqusDnxH0k9tv7duuphq0sUiuiQjyBETRNKXgAeBF9h+WrNa3dm2txnnQ2NASFqVUmR2fmEdSQcCrwduodxT8H3b9zUdV35te4OqAWPKkLTlWMfTYjJqyAhy9KVZlvZI4HGUkUVRBqBmVA02WLa1vaWkuQC2b5XUue4Qkp4DXGH7Dkmvpdyo9znbv68crfNsz6cs17xX7Sx9WB3Yzfbv2jttPyhpYHrnxqTQXnWyPWqnZvsFExsnIjfpRf8+Aexqe1XbM2yvkuJ4id3X3JVtWNgb98G6kUb0JeBOSZsD7wJ+A2TZ4Umk+T389+HFcY/t6yY4UkxhtndsWgvuBPyQcpJ5G3B6sy9iwqVAjn79NS+aD9vngdOAx0k6Avg50MUen+lYMsnZfgC4Pr3Mo2OOA55G+Vt5NLAJpT93xITLFIvo12xJJwPfB+7p7bT9vXqRBovtEyXNAV5IuXT4yo6edCyQdAjwWuC5zZzU5StniqVvdeAaSbOAO3o7be9aL1JMcZva3qS1fZ6ka6uliSktBXL0awZwJ6UlVI+BFMh9krQBcKPt/5X0fODFkm5qeg13STqWTA0rUrqW9Ihyn0FELZdL2s72JQCStgVmV84UU1S6WERMEElXUHrNrkuZZ3c68HTbmWMXE25QejXH5CdpHmXAZXlgI+D3zfY6wC+HjSpHTIiMIEdfJG1IuXnr8bY3lbQZ5aa9j1SONkgetH1/0xHkC7aP7nW06JJ0LJncBq1Xc0wJ6ZoSnZMR5OhLs5Lee4Cv2H5ms+9q25vWTTY4JF0KfA44FNjF9o1d/B5K+j9Kvi7Oj46HaZB6NUdE1JIR5OjXo2zPKovBLXR/rTADah9gf+CIpjhej262T0vHkklswHo1R0RUkQI5+nVLc5NZr4fvq4Cb6kYaLLavBd7W2r6R1k1Rkr5re/ca2YZJx5KIiJjSUiBHvw4AjgE2lvQn4EZKG7BYetavHaCRjiURETGlZQ5yLBFJKwHTbC+onWWyGamrQEREREy8jCDHmCS9c5T9ANj+zIQGimUuHUsiImKqy1LTMZ5VmretgTcDT2re9gcy2rl0afyHTIivAocA9wHYvgr496qJIiIiJlBGkGNMtj8IIOlCYMve1ApJh1MWu4il56DaARrpWBIREVNaCuTo1+OBe1vb9zb7YhytVaIWO0RZgGMzyjtnT2iw0aVjSURETGkpkKNfxwOzJJ3WbL8SOLZenIEyaKtEpWNJRERMaeliEX2TtCXwL83mhbbnto6tbvvWOsliWUjHkoiImKpSIMdSkRZl45O0HXA08DRgBWA6cIftGVWDNUbrWNKTjiURETFVZIpFLC1d6cDQZV+gdIM4ldIV5PXAhlUTDbVK8+9GwDbA6c32LsCsKokiIiIqyAhyLBUZQR6fpNm2t5Z0Ve/GPElzbT+zdra2pmPJy1sdS1YBfmj7uXWTRURETIyMIEdMnDslrQBcIekTlM4QXexFno4lERExpaVAjqUlUyzG9zpKQfzfwDuAtYDdqiYaWTqWRETElJYpFtE3SdMpI4kLT6xs/7459mjb/6iVbRBIOtD2UePt64J0LImIiKksBXL0RdJbgcOAvwIPNrsXLnIR4xtpnnYX5yCPJ/PNIyJisssUi+jXgcBGtv9eO8igkbQXsDewnqTTW4dmAIM46p7pNBERMamlQI5+/QGYXzvEgLqYckPeGsCnW/sXAFdVSfTw5LJTRERMapliEWNqLR7xdEp/3B8C9/SOZ/GIJSPp8ZQewwCzbN9cM89DkSkWEREx2XWxxVR0yyrN2++Bn1JWgFul9RZ9krQHZcGNPYBXA5dKelXdVA9JplhERMSklhHkiAki6Urgxb1RY0mPBc6xvXndZENJ+jBwIXCx7TtGOJ6OJRERMallDnL0RdIZLD73dD4wG/iK7bsnPtXAmTZsSsXf6eZVnBuAvYDPS1oAXERp9fYDgBTHEREx2WUEOfoi6SjgscBJza49gdspRfMM26+rlW1QNKvnbc7Q7+FVtg+ql2p0kp5AmQrybmB125lSExERU0JGkKNfz7a9TWv7DEmX2d5G0jXVUg0WA18Bdmi2jwG2qxdnZJK+BmxC6Xl9EfAq4PKqoSIiIiZQFy/vRjetLGnt3kbz/srN5r11Ig2cF9v+nu13Nm+nAf9aO9QIHgNMB26j9Gm+xfb9dSNFRERMnIwgR7/eBfxc0m8oXQzWA94iaSXguKrJOk7Sm4G3AOtLavc9XgWYWSfV6Gz/G4CkpwEvBc6TNN32k+smi4iImBiZgxx9k/QIYONm8/rcmNcfSasCqwMfAw5uHVrQxRveJO0M/AvwXGA14BLgItvfqBosIiJigqRAjjFJeoHtn0nabaTjtr830Zli2ZL0Bcrc44ts/7l2noiIiImWKRYxnucBPwN2abZ7Z1Rq3k+BPMnY/u/ein+StmRAV/yLiIh4qDKCHH2RtCKwO7Aui06sbPtD1ULFMtGs+Pcp4HzKidC/AO+x/Z2auSIiIiZKRpCjX9+ndDW4HOjNPc7Z1eT0PmCb4Sv+ASmQIyJiSkiBHP16su2X1Q4RE2JQVvyLiIhYJlIgR78ulvQM2/NqB4ll7ieSzmLoin8/qpgnIiJiQmUOcoxJ0jzKVIrlgKcCNwD30NykZ3uzivFiGZG0O/CcZvOiZlGTiIiIKSEFcoxJ0jpjHbf9u4nKEhERETERUiBHBACSFjDyjZe9qwUzJjhSREREFSmQIyIiIiJacmd6RERERERLCuSIiIiIiJYUyBERERERLSmQIyIiIiJa/j+HQ+iKuLH15AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualization of the dataset\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "df = pd.DataFrame(\n",
    "    X,\n",
    "    columns=[\n",
    "        \"age\",\n",
    "        \"sex\",\n",
    "        \"body_mass\",\n",
    "        \"blood_pressure\",\n",
    "        \"total_serum_cholesterol\",\n",
    "        \"low_densisty_lipoproteins\",\n",
    "        \"high_density_lipoproteins\",\n",
    "        \"total_cholesterol\",\n",
    "        \"tryglicerides_level\",\n",
    "        \"blood_sugar_level\",\n",
    "    ],\n",
    ")\n",
    "df[\"target\"] = y\n",
    "corr_mat = df.corr()  # compute the correlation matrix of the features vectors\n",
    "sns.clustermap(corr_mat, vmax=0.8);  # plot a clustered heatmap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Normalization\n",
    "[Normalize](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.normalize.html) the input data using the L2 norm, as well as the target data by the maximum value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO overwrite X,y with normalized versions of the original data\n",
    "X = ...\n",
    "y = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data\n",
    "X_train, X_validation, y_train, y_validation = train_test_split(\n",
    "    X, y, test_size=0.2, random_state=14\n",
    ")\n",
    "X_train, X_test, y_train, y_test = train_test_split(\n",
    "    X_train, y_train, test_size=0.3, random_state=14\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Implementing gradients for Linear Regression\n",
    "\n",
    "The model we consider is $\\hat{y} = w x + c$.\n",
    "\n",
    "We need to calculate the gradients for the training dataset $X$ of the loss\n",
    "$L(w,c) = \\frac{1}{2n}\\sum_{x_i \\in X}(y_i-\\hat{y_i})^2 $\n",
    "\n",
    "What are the gradients?\n",
    "\n",
    "- $\\frac{\\partial L(w,c)}{\\partial w} = ...$\n",
    "- $\\frac{\\partial L(w,c)}{\\partial c} = ...$\n",
    "\n",
    "Implement them below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_w(X, y, w, c):\n",
    "    # outputs parital L(w,c)/ partial w\n",
    "    y_pred = w @ X.T + c\n",
    "    # TODO\n",
    "    return ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_c(X, y, w, c):\n",
    "    # outputs parital L(w,c)/ partial c\n",
    "    y_pred = w @ X.T + c\n",
    "    # TODO\n",
    "    return ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the train and test error.\n",
    "Update the parameters w,c for each gradient step according to the learning rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up stochasic batch gradient descent\n",
    "BATCHSIZE = 16\n",
    "EPOCHS = 400\n",
    "lr = 0.15  # set the learning rate\n",
    "\n",
    "# set initial values for the parameters\n",
    "w = np.random.RandomState(seed=43).randn(X_train.shape[1])\n",
    "c = 0\n",
    "\n",
    "\n",
    "MSE_test = []\n",
    "MSE_train = []\n",
    "\n",
    "# for several epochs (= runs over complete training data)\n",
    "for e in range(EPOCHS):\n",
    "    # TODO\n",
    "    # compute current errors on test and train\n",
    "    MSE_train.append(...)\n",
    "    MSE_test.append(...)\n",
    "\n",
    "    # split the training set into batches\n",
    "    for b in range(0, ((X_train.shape[0] // BATCHSIZE) - 1) * BATCHSIZE, BATCHSIZE):\n",
    "        idxs = np.random.randint(0, X_train.shape[0], BATCHSIZE)\n",
    "\n",
    "        X_batch = X_train[idxs]\n",
    "        y_batch = y_train[idxs]\n",
    "\n",
    "        # TODO\n",
    "        # run a gradient step on every batch\n",
    "        w_ = ...\n",
    "        c_ = ...\n",
    "\n",
    "        w = w_\n",
    "        c = c_\n",
    "\n",
    "# plot the resulting learning curve\n",
    "plt.plot(np.nan_to_num(MSE_train, nan=1000.0), label=\"train\")\n",
    "plt.plot(np.nan_to_num(MSE_test, nan=1000.0), label=\"test\")\n",
    "plt.xlabel(\"epochs\")\n",
    "plt.ylabel(\"MSE\")\n",
    "plt.legend()\n",
    "plt.yscale(\"log\")\n",
    "print(f\"SGD = MSE test {MSE_test[-1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the predictions\n",
    "plt.plot((w @ X_test.T + c), y_test, \".\")\n",
    "plt.xlabel(\"labels predicted by the model\")\n",
    "plt.ylabel(\"true test labels\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Explain in your own words what information the learning curve gives to you."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Loss Curve Quiz\n",
    "\n",
    "What is wrong?\n",
    "![L1](https://developers.google.com/static/machine-learning/testing-debugging/images/metric-curve-ex03.svg)   ![L2](https://developers.google.com/static/machine-learning/testing-debugging/images/metric-curve-ex02.svg)     ![L3](https://developers.google.com/static/machine-learning/testing-debugging/images/metric-curve-ex01.svg\n",
    ")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.4 Comparison with Library solution\n",
    "\n",
    "Why should we (not) use sklearn instead?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the learning curve fot the sklearn Ridge estimator\n",
    "MSE_train = []\n",
    "MSE_test = []\n",
    "for i in range(10, X_train.shape[0]):\n",
    "    lm = linear_model.Ridge(alpha=0.01, fit_intercept=True).fit(\n",
    "        X_train[:i], y_train[:i]\n",
    "    )\n",
    "\n",
    "    # TODO\n",
    "    MSE_test.append(...)\n",
    "    MSE_train.append(...)\n",
    "\n",
    "plt.plot(MSE_train, label=\"train\")\n",
    "plt.plot(MSE_test, label=\"test\")\n",
    "plt.xlabel(\"datapoints\")\n",
    "plt.ylabel(\"MSE\")\n",
    "plt.legend()\n",
    "\n",
    "print(f\"sklearn = MSE test {MSE_test[-1]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lm.predict(X_test), y_test, \".\")\n",
    "plt.xlabel(\"labels predicted by the model\")\n",
    "plt.ylabel(\"true test labels\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. What is the difference between the two methods?\n",
    "2. Should you rather use the Library or your own implementation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7NmeoFdhfQDU"
   },
   "source": [
    "## 2.5 Validation data set\n",
    "1. We did not use the validation data set so far, do we need it at all?\n",
    "2. Are there hyperparamters for SGD?\n",
    "3. How would can you estimate true the generalization error?"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": ".venv",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}