notebooks_en/4_Polynomial_Regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2021 Lorena A. Barba, Tingyu Wang"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial Regression\n",
    "\n",
    "In this fourth lesson of the _Engineering Computations_ module on deep learning, we play with a different model-building method: polynomial regression. \n",
    "We already saw one case in which the observational data could not be fit by a line, and that was the classification problem that we tackled with logistic regression (Lesson 3). \n",
    "Now imagine that your data wiggles about on the $x, y$ plane—where $x$ is the independent variable or feature, and $y$ is the dependent variable. \n",
    "You look at it and think: a curvilinear relationship might work. Good idea. \n",
    "\n",
    "It may surprise you to learn that fitting a polynomial to the data is a special case of multiple linear regression. Yes: _linear_. \n",
    "When we talk about linear models, we mean linear with respect to the parameters!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A special case of multiple linear regression\n",
    "\n",
    "Let's generate some synthetic data using a polynomial function of fourth order, $y = x^4 + x^3 - 4x^2 $, with a bit of added noise.\n",
    "\n",
    "As usuall, we start by loading some needed Python libraries and functions, including the useful tools from `autograd` that you've learned about already."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot\n",
    "from autograd import grad\n",
    "from autograd import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATNUlEQVR4nO3df2xd513H8fcXE+BqA7lV3ZK4LRkoGAZlC7KqoSI01g1XgBZTqWj8UgSVokn8GBKYJUxiAlS1yBICIYSItkEQ+0G0pWnEYF6XMvFDbKu7lKVdalqNrY0TGjOw2MDa0uzLHz4uTuo0uT7n5vg89/2SrHvPc38836Pkfu45z/PceyMzkSSV6evaLkCSNDiGvCQVzJCXpIIZ8pJUMENekgr29W0XsN4NN9yQO3fubLsMSeqUxx577D8yc2yj27ZUyO/cuZP5+fm2y5CkTomIL1zuNodrJKlghrwkFcyQl6SCGfKSVDBDXpIKtqVW10jSsDl6YpHZuQXOLK+wY7THzNQE07vHG3t+Q16SWnL0xCIHjpxk5fwFABaXVzhw5CRAY0HvcI0ktWR2buHFgF+zcv4Cs3MLjfVhyEtSS84sr/TVvhmGvCS1ZMdor6/2zTDkJaklM1MT9LaNXNTW2zbCzNREY3048SpJLVmbXHV1jSQVanr3eKOhfimHaySpYIa8JBXMkJekghnyklQwQ16SCtZIyEfEaER8MCKeiohTEfEDEXF9RDwcEU9Xl9c10Zck6eo1dST/h8BHMvO7gNcAp4D9wPHM3AUcr7YlSddQ7ZCPiG8Bfgh4N0BmfjUzl4E9wKHqboeA6bp9SZL608SR/LcDS8CfRcSJiHhXRLwCuCkzzwJUlzdu9OCI2BcR8xExv7S01EA5kqQ1TYT81wPfD/xJZu4G/oc+hmYy82BmTmbm5NjYWAPlSJLWNBHyp4HTmfnJavuDrIb+8xGxHaC6PNdAX5KkPtQO+cz8d+C5iFj72rQ7gc8Cx4C9Vdte4KG6fUmS+tPUF5T9MvDeiPgG4HPAz7P6BnI4Iu4FngXuaagvSdJVaiTkM/NxYHKDm+5s4vklSZvjJ14lqWCGvCQVzJCXpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQVr5EdDIuLzwJeAC8ALmTkZEdcDfwXsBD4P/GRm/lcT/UmSrk6TR/I/nJmvzcy1X4jaDxzPzF3A8WpbknQNDXK4Zg9wqLp+CJgeYF+SpA00FfIJfDQiHouIfVXbTZl5FqC6vHGjB0bEvoiYj4j5paWlhsqRJEFDY/LAHZl5JiJuBB6OiKeu9oGZeRA4CDA5OZkN1SNJoqEj+cw8U12eAx4Ebgeej4jtANXluSb6kiRdvdohHxGviIhvXrsO/AjwBHAM2FvdbS/wUN2+JEn9aWK45ibgwYhYe773ZeZHIuJR4HBE3As8C9zTQF+SpD7UDvnM/Bzwmg3avwjcWff5JUmb5ydeJalghrwkFcyQl6SCNbVOXpKG0tETi8zOLXBmeYUdoz1mpiaY3j3edlkvMuQlaZOOnljkwJGTrJy/AMDi8goHjpwE2DJB73CNJG3S7NzCiwG/ZuX8BWbnFlqq6KUMeUnapDPLK321t8GQl6RN2jHa66u9DYa8JG3SzNQEvW0jF7X1to0wMzXRUkUv5cSrJG3S2uSqq2skqVDTu8e3VKhfyuEaSSqYIS9JBTPkJalghrwkFayxkI+IkYg4ERF/XW1fHxEPR8TT1eV1TfUlSbo6TR7Jvw04tW57P3A8M3cBx6ttSdI11EjIR8TNwI8B71rXvAc4VF0/BEw30Zck6eo1dST/B8BvAF9b13ZTZp4FqC5vbKgvSdJVqh3yEfHjwLnMfGyTj98XEfMRMb+0tFS3HEnSOk0cyd8BvDkiPg98AHhDRPwl8HxEbAeoLs9t9ODMPJiZk5k5OTY21kA5kqQ1tUM+Mw9k5s2ZuRN4C/BIZv4scAzYW91tL/BQ3b4kSf0Z5Dr5B4A3RcTTwJuqbUnSNdToF5Rl5seBj1fXvwjc2eTzS5L64ydeJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekgjX63TVtOXpikdm5Bc4sr7BjtMfM1ATTu8fbLkuSWtf5kD96YpEDR06ycv4CAIvLKxw4chLAoJc09Do/XDM7t/BiwK9ZOX+B2bmFliqSpK2j8yF/Znmlr3ZJGiadD/kdo72+2iVpmDTxQ97fFBGfioh/iYgnI+K3q/brI+LhiHi6uryufrkvNTM1QW/byEVtvW0jzExNDKI7SeqUJo7kvwK8ITNfA7wWuCsiXgfsB45n5i7geLXduOnd49x/922Mj/YIYHy0x/133+akqyTRwOqazEzgy9XmtuovgT3A66v2Q6z+LODb6/a3kend44a6JG2gkTH5iBiJiMeBc8DDmflJ4KbMPAtQXd54mcfui4j5iJhfWlpqohxJUqWRkM/MC5n5WuBm4PaI+N4+HnswMyczc3JsbKyJciRJlUZX12TmMqvDMncBz0fEdoDq8lyTfUmSrqyJ1TVjETFaXe8BbwSeAo4Be6u77QUeqtuXJKk/TXytwXbgUESMsPqmcTgz/zoi/hk4HBH3As8C9zTQlySpD02srvkMsHuD9i8Cd9Z9fknS5nX+E6+SpMsz5CWpYIa8JBXMkJekghnyklQwQ16SCmbIS1LBDHlJKpghL0kFM+QlqWCGvCQVzJCXpII18S2UktSaoycWmZ1b4MzyCjtGe8xMTfhzoOsY8pI66+iJRQ4cOcnK+QsALC6vcODISQCDvmLIS+qs2bmFFwN+zcr5C8zOLVx1yJd+JmDIS+qsM8srfbVfahjOBJr4+b9bIuLvIuJURDwZEW+r2q+PiIcj4unq8rr65UrS/9sx2uur/VIvdyZQiiZW17wA/FpmfjfwOuAXI+LVwH7geGbuAo5X25LUmJmpCXrbRi5q620bYWZq4qoeX/dMoAtqh3xmns3MT1fXvwScAsaBPcCh6m6HgOm6fUnSetO7x7n/7tsYH+0RwPhoj/vvvu2qh1rqngl0QaNj8hGxk9Xfe/0kcFNmnoXVN4KIuPEyj9kH7AO49dZbmyxH0hCY3j2+6fHzmamJi8bkob8zgS5o7MNQEfFK4EPAr2bmf1/t4zLzYGZOZubk2NhYU+VI0hXVPRPogkaO5CNiG6sB/97MPFI1Px8R26uj+O3AuSb6kqQm1TkT6IImVtcE8G7gVGb+/rqbjgF7q+t7gYfq9iVJ6k8TR/J3AD8HnIyIx6u23wQeAA5HxL3As8A9DfQlSepD7ZDPzH8E4jI331n3+SVJm+e3UEpSwQx5SSqYIS9JBTPkJalghrwkFcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQUz5CWpYIa8JBXMkJekgjX6G6+S+nf0xCKzcwucWV5hx2iPmamJon+pSNeWIb8F+CIfXkdPLF70Q9KLyyscOHISwP8DakQjwzUR8Z6IOBcRT6xruz4iHo6Ip6vL65roqzRrL/LF5RWS/3+RHz2x2HZpugZm5xZeDPg1K+cvMDu30FJFKk1TY/J/Dtx1Sdt+4Hhm7gKOV9u6hC/y4XZmeaWvdqlfjYR8Zv498J+XNO8BDlXXDwHTTfRVGl/kw23HaK+vdqlfg1xdc1NmngWoLm/c6E4RsS8i5iNifmlpaYDlbE2+yIfbzNQEvW0jF7X1to0wMzXRUkUqTetLKDPzYGZOZubk2NhY2+Vcc77Ih9v07nHuv/s2xkd7BDA+2uP+u29z0lWNGeTqmucjYntmno2I7cC5AfbVWWsvZlfXDK/p3eP+e2tgBhnyx4C9wAPV5UMD7KvTfJFLGpRGQj4i3g+8HrghIk4D72Q13A9HxL3As8A9TfQ1CK5Tl1SqRkI+M3/qMjfd2cTzD1IJH0bxTUrS5bQ+8dq2rq9T98NUkl7O0Id819epd/1NStJgDX3Id32detffpCQN1tCHfNfXqXf9TUo6emKROx54hFft/zB3PPCIQ40NG/qQ7/qHUbr+JqX6uhySzikNnl81TLfXqfthquHW9dVhLzen1IX6u8CQL0CX36RUT9dD0jmlwTPkG+A6dbWl6yG5Y7TH4ga1OqfUHEO+pq6fLjeh629yXa6/6yE5MzVx0esHnFNq2tBPvNY17OvUuz5x1vX6uz7x3vWFD13gkXxNXT9drqvrY8Jdr38rTLzXPRNyTmmwDPmaun66XFfX3+S6Xj+0G5IOV259DtfU1PXT5bq6/mGsrtfftmEfruwCQ76mYR9T7PqbXBP1d/nDSHWVcCZUOodrGjDMY4pbYUy4jrr1D/twxbAPV3ZBZOZgO4i4C/hDYAR4V2Y+cLn7Tk5O5vz8/EDrkZp0xwOPbBhy46M9/mn/G1qo6Nq69E0OVs+EhulsdiuIiMcyc3Kj2wZ6JB8RI8AfA28CTgOPRsSxzPzsIPuVrpVhH67o+pncMBj0cM3twDOZ+TmAiPgAsAcw5FUEhyuGe7iyCwY98ToOPLdu+3TVJhWh6xPPMNwTx8Ng0EfysUHbRZMAEbEP2Adw6623DrgcqVldH64Y9onjYTDokD8N3LJu+2bgzPo7ZOZB4CCsTrwOuB4VqO3vnunycEXXP/GrKxt0yD8K7IqIVwGLwFuAnx5wnxoiHonWM+wTx8NgoGPymfkC8EvAHHAKOJyZTw6yTw0XP3FZj5/4Ld/AP/GamX+Tmd+Zmd+RmfcNuj8NF49E6ylh4lgvz681UKd5JFrPsH8txzDwaw3UujoTp/7oRH1dnjjWlRnyalXdidOuL2GUBs2QV6uaWMLnkah0eY7Jq1VOnEqDZcirVU6cSoNlyKtVLuGTBssxebXKiVNpsAx5tc6JU2lwHK6RpIIZ8pJUMENekgpmyEtSwQx5SSqYIS9JBTPkJalgtUI+Iu6JiCcj4msRMXnJbQci4pmIWIiIqXplSpI2o+6HoZ4A7gb+dH1jRLya1d9z/R5gB/CxiPjOzLzw0qeQJA1KrSP5zDyVmRv9mOYe4AOZ+ZXM/DfgGeD2On1Jkvo3qDH5ceC5ddunq7aXiIh9ETEfEfNLS0sDKkeShtMVh2si4mPAt25w0zsy86HLPWyDttzojpl5EDgIMDk5ueF9JEmbc8WQz8w3buJ5TwO3rNu+GTizieeRJNUwqOGaY8BbIuIbI+JVwC7gUwPqS5J0GXWXUP5ERJwGfgD4cETMAWTmk8Bh4LPAR4BfdGWNJF17tZZQZuaDwIOXue0+4L46zy9JqsdPvEpSwQx5SSqYIS9JBfM3XsXRE4v+kLZUKEN+yB09sciBIydZOb+6+GlxeYUDR04CGPRSARyuGXKzcwsvBvyalfMXmJ3b6CuJJHWNIT/kziyv9NUuqVsM+SG3Y7TXV7ukbjHkh9zM1AS9bSMXtfW2jTAzNdFSRZKa5MTrkFubXHV1jVQmQ15M7x431KVCOVwjSQUz5CWpYIa8JBXMkJekghnyklSwyNw6v50dEUvAF2o8xQ3AfzRUTptK2Q9wX7aiUvYD3Jc135aZYxvdsKVCvq6ImM/MybbrqKuU/QD3ZSsqZT/AfbkaDtdIUsEMeUkqWGkhf7DtAhpSyn6A+7IVlbIf4L5cUVFj8pKki5V2JC9JWseQl6SCFRXyEfG7EfGZiHg8Ij4aETvarmmzImI2Ip6q9ufBiBhtu6bNioh7IuLJiPhaRHRuuVtE3BURCxHxTETsb7uezYqI90TEuYh4ou1a6oqIWyLi7yLiVPV/621t17QZEfFNEfGpiPiXaj9+u/E+ShqTj4hvycz/rq7/CvDqzHxry2VtSkT8CPBIZr4QEb8HkJlvb7msTYmI7wa+Bvwp8OuZOd9ySVctIkaAfwXeBJwGHgV+KjM/22phmxARPwR8GfiLzPzetuupIyK2A9sz89MR8c3AY8B01/5dIiKAV2TmlyNiG/CPwNsy8xNN9VHUkfxawFdeAXT2HSwzP5qZL1SbnwBubrOeOjLzVGZ29ZfBbweeyczPZeZXgQ8Ae1quaVMy8++B/2y7jiZk5tnM/HR1/UvAKaBzP4qQq75cbW6r/hrNraJCHiAi7ouI54CfAX6r7Xoa8gvA37ZdxJAaB55bt32aDoZJySJiJ7Ab+GTLpWxKRIxExOPAOeDhzGx0PzoX8hHxsYh4YoO/PQCZ+Y7MvAV4L/BL7Vb78q60L9V93gG8wOr+bFlXsy8dFRu0dfYMsTQR8UrgQ8CvXnIm3xmZeSEzX8vq2frtEdHoUFrnfv4vM994lXd9H/Bh4J0DLKeWK+1LROwFfhy4M7f45Ekf/y5dcxq4Zd32zcCZlmrROtUY9oeA92bmkbbrqSszlyPi48BdQGOT4507kn85EbFr3eabgafaqqWuiLgLeDvw5sz837brGWKPArsi4lUR8Q3AW4BjLdc09KoJy3cDpzLz99uuZ7MiYmxt5VxE9IA30nBulba65kPABKsrOb4AvDUzF9utanMi4hngG4EvVk2f6PBKoZ8A/ggYA5aBxzNzqtWi+hARPwr8ATACvCcz72u3os2JiPcDr2f1K22fB96Zme9utahNiogfBP4BOMnq6x3gNzPzb9qrqn8R8X3AIVb/b30dcDgzf6fRPkoKeUnSxYoarpEkXcyQl6SCGfKSVDBDXpIKZshLUsEMeUkqmCEvSQX7PxmPnhXp37F8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "numpy.random.seed(0)    # fix seed for reproducibility\n",
    "x = numpy.linspace(-3, 3, 20)\n",
    "y = x**4 + x**3 - 4*x**2 + 8*numpy.random.normal(size=len(x))\n",
    "pyplot.scatter(x, y);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we were only given the data (not the function that generated them), and our goal is to fit a curve to these data. The nonlinear relationship between $x$ and $y$ suggests that a linear regression will fail. Intuitively, using polynomial functions may first come to your mind.\n",
    "\n",
    "Let's write the model as a $d$th-order polynomial on $x$, the only feature:\n",
    "\n",
    "$$\n",
    "\\hat{y} = w_0 + w_1 x + w_2 x^2 + \\cdots + w_d x^d, \n",
    "$$\n",
    "\n",
    "where $w$ denotes the weights. Keep in mind that in the model-fitting context, the objective is always to find the optimal values of these weights given $x$ and $y$. When viewed from a different perspective, the model above is just a linear combination of the weights. In fact, by creating polynomial features of $x$, namely, letting $x_i = x^i$, the model becomes:\n",
    "\n",
    "$$\n",
    "\\hat{y} = w_0 + w_1 x_1 + w_2 x_2 + \\ldots + w_d x_d.\n",
    "$$\n",
    "\n",
    "As you can see, the polynomial regression model is identical to multiple linear regression, with the matrix form being also $\\hat{\\mathbf{y}} = X\\mathbf{w}$, and the only gap is forming the matrix $X$ using the powers of $x$.\n",
    "\n",
    "Suppose we want to fit our data with a 3rd degree polynomial function; let's write a function to create these polynomial features first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(20, 4)\n"
     ]
    }
   ],
   "source": [
    "degree = 3\n",
    "\n",
    "def polynomial_features(x, degree):\n",
    "    \"\"\" Generate polynomial features for x.\"\"\"\n",
    "    \n",
    "    X = numpy.empty((len(x), degree+1))\n",
    "    for i in range(degree+1):\n",
    "        X[:,i] = x**i\n",
    "    return X\n",
    "\n",
    "X = polynomial_features(x, degree)\n",
    "print(X.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Note\n",
    "\n",
    "> Unsurprisingly, **scikit-learn** offers a counterpart: [`PolynomialFeatures()`](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PolynomialFeatures.html).\n",
    "When your original data come with multiple features (e.g., $x_1$ and $x_2$), the polynomial features for regression will involve interaction terms (e.g., $x_1 x_2$ for 2nd-order polynomials).\n",
    "In that case, it is handy to use this function to generate all terms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scale the data, train the model\n",
    "\n",
    "Recall that for multiple linear regression, we should normalize each feature to the same scale.\n",
    "As in the previous lesson, let's use `MinMaxScaler()` to scale all features to $[0,1]$, except for the first column:\n",
    "$x_0$ is set to 1 for all entries, since $w_0$ represents the intercept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "\n",
    "min_max_scaler = MinMaxScaler()\n",
    "X_scaled = min_max_scaler.fit_transform(X)\n",
    "X_scaled[:,0] = 1   # the column for intercept"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When scaling the matrix of polynomial features $X$ above, we used:\n",
    "```python\n",
    "X_scaled = min_max_scaler.fit_transform(X)\n",
    "```\n",
    "The function `fit_transform()` is actually a combination of two steps:\n",
    "- [`fit()`](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.MinMaxScaler.html#sklearn.preprocessing.MinMaxScaler.fit): compute $x_{\\min}$ and $x_{\\max}$ values for each feature and save the information to the variable `min_max_scaler`. You can access them via `min_max_scaler.data_min_` and `min_max_scaler.data_max_`.\n",
    "- [`transform()`](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.MinMaxScaler.html#sklearn.preprocessing.MinMaxScaler.transform) use `min_max_scaler` to scale `X`.\n",
    "\n",
    "This will be helpful to remember a bit later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can reuse the same model and loss function from Lesson 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_regression(params, X):\n",
    "    '''\n",
    "    The linear regression model in matrix form.\n",
    "    Arguments:\n",
    "      params: 1D array of weights for the linear model\n",
    "      X     : 2D array of input values\n",
    "    Returns:\n",
    "      1D array of predicted values\n",
    "    '''\n",
    "    return numpy.dot(X, params)\n",
    "\n",
    "def mse_loss(params, model, X, y):\n",
    "    '''\n",
    "    The mean squared error loss function.\n",
    "    Arguments:\n",
    "      params: 1D array of weights for the linear model\n",
    "      model : function for the linear regression model\n",
    "      X     : 2D array of input values\n",
    "      y     : 1D array of predicted values\n",
    "    Returns:\n",
    "      float, mean squared error\n",
    "    '''\n",
    "    y_pred = model(params, X)\n",
    "    return numpy.mean( numpy.sum((y-y_pred)**2) )\n",
    "\n",
    "gradient = grad(mse_loss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember, \"training\" a model simply means finding the best parameters by minimizing a loss function. \n",
    "We'll choose both a maximum number of iterations in the optimization loop, and an exit criterion based on the norm of the gradient at the current iteration. \n",
    "(If this is quite small, when multiplied by the also-small learning rate, the parameters will change very little.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration    0, loss = 5434.768, mae = 11.057\n",
      "iteration  100, loss = 1300.477, mae = 6.885\n",
      "iteration  200, loss = 1281.309, mae = 6.864\n",
      "iteration  300, loss = 1272.990, mae = 6.808\n",
      "iteration  400, loss = 1267.448, mae = 6.760\n",
      "iteration  500, loss = 1263.750, mae = 6.722\n",
      "iteration  600, loss = 1261.282, mae = 6.690\n",
      "iteration  700, loss = 1259.636, mae = 6.664\n",
      "iteration  800, loss = 1258.537, mae = 6.643\n",
      "iteration  900, loss = 1257.804, mae = 6.625\n",
      "iteration 1000, loss = 1257.314, mae = 6.611\n",
      "iteration 1100, loss = 1256.988, mae = 6.600\n",
      "iteration 1200, loss = 1256.770, mae = 6.590\n",
      "iteration 1300, loss = 1256.625, mae = 6.583\n",
      "iteration 1400, loss = 1256.528, mae = 6.576\n",
      "iteration 1500, loss = 1256.463, mae = 6.571\n",
      "iteration 1600, loss = 1256.420, mae = 6.567\n",
      "iteration 1700, loss = 1256.391, mae = 6.564\n",
      "iteration 1800, loss = 1256.372, mae = 6.561\n",
      "iteration 1900, loss = 1256.359, mae = 6.558\n",
      "iteration 2000, loss = 1256.350, mae = 6.557\n",
      "iteration 2100, loss = 1256.345, mae = 6.555\n",
      "iteration 2200, loss = 1256.341, mae = 6.554\n",
      "iteration 2300, loss = 1256.338, mae = 6.553\n",
      "iteration 2400, loss = 1256.337, mae = 6.552\n",
      "iteration 2500, loss = 1256.336, mae = 6.551\n",
      "iteration 2600, loss = 1256.335, mae = 6.551\n",
      "iteration 2700, loss = 1256.334, mae = 6.550\n",
      "iteration 2800, loss = 1256.334, mae = 6.550\n",
      "iteration 2900, loss = 1256.334, mae = 6.550\n"
     ]
    }
   ],
   "source": [
    "max_iter = 3000\n",
    "alpha = 0.01\n",
    "params = numpy.zeros(X_scaled.shape[1])\n",
    "descent = numpy.ones(X_scaled.shape[1])\n",
    "i = 0\n",
    "\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "while numpy.linalg.norm(descent) > 0.01 and i < max_iter:\n",
    "    descent = gradient(params, linear_regression, X_scaled, y)\n",
    "    params = params - descent * alpha\n",
    "    loss = mse_loss(params, linear_regression, X_scaled, y)\n",
    "    mae = mean_absolute_error(y, X_scaled@params)\n",
    "    if i%100 == 0:\n",
    "        print(f\"iteration {i:4}, {loss = :.3f}, {mae = :.3f}\")\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print out the weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-22.51572398,   6.75930601,  41.30788709,  30.0105898 ])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first term is the intercept, and the rest are the weights of each **scaled** polynomial feature.\n",
    "\n",
    "Although our model has multiple weights, there's only one feature in our original data. \n",
    "Therefore, if we are given new values of $x$ to predict $y$, we need to first create the polynomial features, scale them to $[0,1]$ using the **same scaling function**, and then multiply by the weights. Since we have used the min-max scaling, it is important to use the same $x_{\\min}$ and $x_{\\max}$ of the training data to scale new data.\n",
    "\n",
    "Recall the explanation above for the scaling step. We now only need to call `min_max_scaler.transform()` to scale new data. Such design occurs very often in scikit-learn.\n",
    "\n",
    "With that in mind, let's plot the fitted curve together with the data. We generate some new values of $x$, `xgrid`, and predict $y$ at these locations. Don't forget to repeat the procedure of creating polynomial features and scaling them. Ponder over this code for a moment to wrap your head around it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAqzElEQVR4nO3deXxU1f3/8dcnC4QEwhIChIQsLLKvCTsqCihWCm64hV+pG1qrdrEqFqu2lqqV+tXaKqaiYolYNyqiRZBFQNawyS5LFsKWEAgkhOzn98dMMGACSeZO7szk83w8eNzcOzP3fiaa95w599x7xBiDUkop3+RndwFKKaXcR0NeKaV8mIa8Ukr5MA15pZTyYRrySinlwwLsLqCy1q1bm9jYWLvLUEopr7Jx48bjxpjwqh7zqJCPjY0lJSXF7jKUUsqriEh6dY9pd41SSvkwDXmllPJhGvJKKeXDPKpPviolJSVkZmZSWFhodyleLygoiKioKAIDA+0uRSlVTzw+5DMzM2nWrBmxsbGIiN3leC1jDDk5OWRmZhIXF2d3OUqpeuLx3TWFhYWEhYVpwLtIRAgLC9NvREp5muRkiI0FPz/HMjnZ0t17fEse0IC3iP4elfIwyckwZQoUFDjW09Md6wCJiZYcwuNb8kop5bOmTfsh4CsUFDi2W0RDvp4tX76ccePGATB//nxeeOGFap+bm5vL66+/XutjPPvss8yYMaPONSql6klGRu2214GGvEXKyspq/Zrx48czderUah+va8grpbxEdHTttteBhnwNpKWl0a1bNyZPnkyfPn245ZZbKCgoIDY2lj/96U+MGDGCjz76iEWLFjF06FAGDBjAxIkTyc/PB2DhwoV069aNESNG8Omnn57b77vvvstDDz0EwLFjx7jxxhvp27cvffv2ZfXq1UydOpX9+/fTr18/HnvsMQBeeuklBg4cSJ8+fXjmmWfO7Wv69Ol07dqV0aNHs2fPnnr87Sil6mz6dAgOPn9bcLBju0W84sRrhT9+voOdh09bus8e7UN55qc9L/m8PXv2MGvWLIYPH87dd999roUdFBTEqlWrOH78ODfddBNff/01ISEhvPjii7z88ss8/vjj3HfffSxdupTOnTtz2223Vbn/Rx55hCuvvJJ58+ZRVlZGfn4+L7zwAtu3b2fLli0ALFq0iL1797J+/XqMMYwfP54VK1YQEhLCBx98wObNmyktLWXAgAHEx8db9jtSSrlJxcnVadMcXTTR0Y6At+ikK3hZyNupQ4cODB8+HIBJkybx97//HeBcaK9du5adO3eee05xcTFDhw5l9+7dxMXF0aVLl3OvTUpK+tH+ly5dynvvvQeAv78/zZs35+TJk+c9Z9GiRSxatIj+/fsDkJ+fz969e8nLy+PGG28k2NkiGD9+vNVvXynlLomJlob6hbwq5GvS4naXC4cfVqyHhIQAjouNxowZw9y5c8973pYtWywbumiM4cknn+T+++8/b/srr7yiwyOVUlXSPvkaysjIYM2aNQDMnTuXESNGnPf4kCFD+Pbbb9m3bx8ABQUFfP/993Tr1o3U1FT2799/7rVVGTVqFG+88QbgOIl7+vRpmjVrRl5e3rnnXHvttbz99tvn+voPHTpEVlYWV1xxBfPmzePs2bPk5eXx+eefW/vmlVJeS0O+hrp3787s2bPp06cPJ06c4Be/+MV5j4eHh/Puu+9yxx130KdPH4YMGcLu3bsJCgoiKSmJ66+/nhEjRhATE1Pl/l999VWWLVtG7969iY+PZ8eOHYSFhTF8+HB69erFY489xjXXXMOdd97J0KFD6d27N7fccgt5eXkMGDCA2267jX79+nHzzTdz+eWX18evRCnlBcQYY3cN5yQkJJgLJw3ZtWsX3bt3t6kih7S0NMaNG8f27dttrcMKnvD7VEpZS0Q2GmMSqnpMW/JKKeXDLAl5EWkhIh+LyG4R2SUiQ0WklYgsFpG9zmVLK45lh9jYWJ9oxSulGh6rWvKvAguNMd2AvsAuYCqwxBjTBVjiXFdKKVWPXA55EQkFrgBmARhjio0xucAEYLbzabOBG1w9llJKqdqxoiXfEcgG3hGRzSLyloiEAG2NMUcAnMs2Vb1YRKaISIqIpGRnZ1tQjlJKqQpWhHwAMAB4wxjTHzhDLbpmjDFJxpgEY0xCeHi4BeUopZSqYEXIZwKZxph1zvWPcYT+MRGJAHAusyw4li30bpBKKW/lcsgbY44CB0Wkq3PTKGAnMB+Y7Nw2GfjM1WPZpbqQr8vthZVSqj5ZNbrmYSBZRL4D+gF/AV4AxojIXmCMc9393DBfYuVb/g4cOJCrrrqKO++8k969e5OWlkavXr3OPXfGjBk8++yzAOzfv5+xY8cSHx/P5Zdfzu7du12uRSmlasOSG5QZY7YAVV1tNcqK/deYm+ZLrHzL3+XLl3P99dezfft24uLiSEtLq/Z1U6ZMYebMmXTp0oV169bx4IMPsnTp0jrXoZRSteVVd6G8pIvNl2jhrTwHDRpEXFzcRZ+Tn5/P6tWrmThx4rltRUVFltWglFI14VshXw/zJcIPtxcGCAgIoLy8/Nx6YWEhAOXl5bRo0eLchB9KKWUH37p3jZvmS7zwlr+VtW3blqysLHJycigqKmLBggUAhIaGEhcXx0cffQQ47gW/detWl+pQSqna8q2Qd9N8iRfe8reywMBAnn76aQYPHsy4cePo1q3buceSk5OZNWsWffv2pWfPnnz2mdcOMFJKeSnfu9VwcrJb50v0dnqrYaU8z6MfbqVjeAi/vKpznV7fsG41nJgIaWlQXu5YasArpTzY7qOn+WRTJkWl5Zd+ch34XsgrpZQXeW3pPkIa+XP38Fi37N8rQt6TupS8mf4elfIs+7Ly+HLbESYPi6VFcCO3HMPjQz4oKIicnBwNKBcZY8jJySEoKMjuUpRSTv9Yuo+gAH/uGXHx625c4fHj5KOiosjMzERvQ+y6oKAgoqKi7C5DKQWkHj/D/K2HuffyjoQ1bey243h8yAcGBl7y6lKllPI2/1y2j0B/P+67vKNbj+Px3TVKKeVrDp4oYN7mQ9w5OJrwZu5rxYOGvFJK1bvXl+/D30944MpObj+WhrxSStWjQ7ln+XhjJrcldKBtqPsHQmjIK6VUPZq5fD8AD4x0fyseNOSVUqreHD1VyH82HOSW+CgiWzSpl2NaMrpGRNKAPKAMKDXGJIhIK+A/QCyQBtxqjDlpxfGUUsobvbliP2XG8ODIut2jpi6sbMlfZYzpV+kmOVOBJcaYLsAS57pSSjVIWXmFvL8ugxv7R9KhVfClX2ARd3bXTABmO3+eDdzgxmMppZRHe2tlKiVl5XW+02RdWRXyBlgkIhtFxDmpKm2NMUcAnMs2Vb1QRKaISIqIpOhVrUopX5STX8S/16QzoV8kca1DLv0CC1l1xetwY8xhEWkDLBaR3TV9oTEmCUgCx/3kLapHKaU8xqxVqRSWltV7Kx4saskbYw47l1nAPGAQcExEIgCcyywrjqWUUt4kt6CY99akc33vCDq3aVrvx3c55EUkRESaVfwMXANsB+YDk51Pmwzo3HdKqQbn7W/TyC8q5aGr678VD9Z017QF5olIxf7eN8YsFJENwIcicg+QAUy04FhKKeU1TheW8M63qYzt2Y5u7UJtqcHlkDfGHAD6VrE9Bxjl6v6VUspbzf42jbxC+1rxoFe8KqWUW+QXlTLr21RGd29Dr8jmttWhIa+UUm7wzqpUcgtKePjqLrbWoSGvlFIWy8kv4s0VB7i2Z1v6dmhhay0a8kop5YrkZIiNBT8/xzI5mX8s28fZkjIeH9vN7uo05JVSqs6Sk2HKFEhPB2MgPZ2M3z3FnNWp3JrQgU7h9T8u/kIa8kopVVfTpkFBwXmb/pZwC/4lJfx6tL198RU05JVSqq4yMs5b3d6mI5/1HMk9G/5bL7M+1YSGvFJK1VV09HmrL478OS0LTnH/4XU2FfRjGvJKKVVX06dDsOPe8Kti+rIybgAPpcwj9Nk/2FzYD6y6C6VSSjU8iYkAlE97ihdG3kXkmRwm/eKGc9s9gYa8Ukq5IjGRBb2uYvvczfzfbX1p3D/K7orOo901SinlguLScmZ8tYfuEaFM6Btpdzk/oiGvlFIumLs+g4wTBUy9rht+fmJ3OT+iIa+UUnWUV1jC35fsZVinMK7o0trucqqkIa+UUnX0r5Wp5JwpZup13XDOqeFxLAt5EfEXkc0issC53kpEFovIXueypVXHUkopu2XlFfLWygOM6xNBn6gWdpdTLStb8r8CdlVanwosMcZ0AZY415VSyie8tmQfxaXl/O6arnaXclGWhLyIRAHXA29V2jwBmO38eTZwgxXHUkopu6UeP8Pc9RncOTia2NYhdpdzUVa15F8BHgfKK21ra4w5AuBctrHoWEopZasZX+2hcYCf7ROC1ITLIS8i44AsY8zGOr5+ioikiEhKdna2q+UopZRbbTmYyxfbjnDfFR0Jb9bY7nIuyYqW/HBgvIikAR8AV4vIHOCYiEQAOJdZVb3YGJNkjEkwxiSEh4dbUI5SSrmHMYbnv9xF66aNuPfyjnaXUyMuh7wx5kljTJQxJha4HVhqjJkEzAcmO582GfjM1WMppZSdPv/uCOtST/CbMZfRtLF33BXGnePkXwDGiMheYIxzXSmlvFJ+USnTv9hJ78jm3D4w+tIv8BCWfhQZY5YDy50/5wCjrNy/UkrZ5bUlezl2uoiZk+Lx98DbF1RHr3hVSqlL2JeVx6xVqdyW0IH+0d51XaeGvFJKXYQxhqc/20FwI38eH+vZFz5VRUNeKaUu4ottR1i9P4fHru1KWFPPHzJ5IQ15pZSqxpmiUv68YBc924dy5+AYu8upE+8YA6SUUjZ4bek+jp4u5J+JA7zqZGtl2pJXSqkq7M/OZ9aqA0yMjyI+xrtOtlamIa+UUhcwxvDs/B00CfTnieu62V2OSzTklVLqAv/bfpSVe4/z6DVdae2FJ1sr05BXSqlKCopL+fOCnXSPCCVxsPdc2VodDXmllKrkH0v3cfhUIc9N6EmAv/dHpPe/A6WUssiB7Hz+tfIANw2IJCG2ld3lWMI3Qj45GWJjwc/PsUxOtrsipZSXMcbwzPwdBAX48+R13e0uxzLeP04+ORmmTIGCAsd6erpjHSAx0b66lFJe5asdjpOtT4/r4RWTgdSU97fkp037IeArFBQ4tiulVA2cLS7juQW76NauGT8b6p1XtlbH+1vyGRm1266UUhf426I9HMo9y3+mDPGJk62Vef+7ia5miFN125VSqpKN6SeY9W0qiYOjGdwxzO5yLGfFRN5BIrJeRLaKyA4R+aNzeysRWSwie51L91wXPH06BAefvy042LFdKaUuorCkjMc++o72zZvw5E9852RrZVa05IuAq40xfYF+wFgRGQJMBZYYY7oAS5zr1ktMhKQkiIkBEccyKUlPuiqlLmnGV3s4cPwMf72lj9fM2VpbLr8rY4wB8p2rgc5/BpgAjHRun41jWsAnXD1elRITNdSVUrWSkvZDN83wzq3tLsdtLOmTFxF/EdkCZAGLjTHrgLbGmCMAzmWbal47RURSRCQlOzvbinKUUuqiCkvKeOxj3+6mqWBJyBtjyowx/YAoYJCI9KrFa5OMMQnGmITw8HArylFKqYua8dUeUo+f4SUf7qapYOnoGmNMLo5umbHAMRGJAHAus6w8llJK1UXlbpphPtxNU8GK0TXhItLC+XMTYDSwG5gPTHY+bTLwmavHUkopV5wtbjjdNBWs+J4SAcwWEX8cHxofGmMWiMga4EMRuQfIACZacCyllKqzGYsc3TTv3zvY57tpKlgxuuY7oH8V23OAUa7uXymlrJCSdoK3v01l0pCG0U1TwfuveFVKqUuo3E0z1YfuMFkTDeP7ilKqQWuI3TQVtCWvlPJpGxpoN00FDXmllM86W1zG4x9/R2SLJj41EUhtNKzvLUqpBuWlr37opglpYN00FbQlr5TySUt2HePtb1P52dCYBtlNU0FDXinlcw7lnuXRj7bSPSKU3zeQi56qoyGvlPIpJWXlPPz+JkpKy3k9cQBBgf52l2SrhtlJpZTyWTO+2sOmjFz+fkd/4lqH2F2O7bQlr5TyGUt2HePNFQdIHBzN+L7t7S7HI2jIK6V8QkU/fI+IUP4wrofd5XgMDXmllNer6IcvLTP8U/vhz6Mhr5TybsnJzLjpt2zKyOX5pW8S99V/7a7Io2jIK6W8V3IyS/4ykze7j2HSpi/46ap5MGUKJCfbXZnH0JBXSnmtQ9Nn8OjoB+lxbD9PLX3LsbGgAKZNq/lOkpMhNhb8/BxLH/uA0CGUSimvVFJWzsPxiZT6BfDPz14kqKzkhwczMmq2k+RkR8u/oMCxnp7uWAdITLS2YJtYMf1fBxFZJiK7RGSHiPzKub2ViCwWkb3OZUvXy1VKKYcZX+1hU2R3nl/4GnEnD5//YHR0zXYybdoPAV+htt8EPJwV3TWlwKPGmO7AEOCXItIDmAosMcZ0AZY415VSymUV4+EntSripxkbz38wOBimT6/Zjqpr8df0m4AXcDnkjTFHjDGbnD/nAbuASGACMNv5tNnADa4eSymlMk8WnBsP/9RvJkBSEsTEgIhjmZRU866W6lr8Nf0m4AUsPfEqIrE45ntdB7Q1xhwBxwcB0Kaa10wRkRQRScnOzrayHKWUj8krLOGed1MoK680Hj4xEdLSoLzcsaxNX/r06Y6Wf2W1+SbgBSwLeRFpCnwC/NoYc7qmrzPGJBljEowxCeHh4VaVo5TyMaVl5Tz0/mb2Z+czc1K8NfelSUx07ZuAF7Ak5EUkEEfAJxtjPnVuPiYiEc7HI4AsK45VlZz8Im5+YzXL97jtEEopGxlj+OPnO/nm+2yeu6EXw628P7wr3wS8gBWjawSYBewyxrxc6aH5wGTnz5OBz1w9VnWOnCokJ7+In7+zgZ+9vZ49R/PcdSillA3eXZ3Gv9emc/8VHbljkO/0l9cHMca4tgOREcBKYBtQ7tz8exz98h8C0UAGMNEYc+Ji+0pISDApKSl1qqO4tJx/r03n1a+/J7+olNsHRfOb0ZcR3qxxnfanlPIMS3Yd4773UhjdvS0zJ8Xj5yd2l+RxRGSjMSahysdcDXkruRLyFU6eKebVJXuZszadoEB/HryqE3cPj9MbFinlhXYePs0tM1fTKbwp/7l/CMGN9PrNqlws5H3utgYtQxrx7PiefPWbKxjSMYy/LtzDqL99w+dbD+NJH2hKqYs7drqQe2ZvoHmTQN6anKABX0c+F/IVOoU35a3JCbx/72BCmwTy8NzN3PTGajZlnLS7NKXUJRQUl3Lv7BROny1h1uSBtA0Nsrskr+WzIV9hWOfWLHh4BH+9uQ+ZJ89y0+ureXjuZjJPFlz6xUqpeldebvj1B1vYcfgUr93Znx7tQ+0uyav5fMgD+PsJtw7swPLfjeThqzuzeOdRrv7bN/x14W7yi0rtLk8pVcmLC3ezaOcxnh7Xg6u7tbW7HK/XIEK+QkjjAB69pitLHx3J9b0jeH35fka+tIy56zMoK9f+eqXsNnd9Bm+uOMDkoTH8fHic3eX4hAYV8hXat2jC/93Wj//+cjixYSE8+ek2rv/7SlbtPW53aUo1WKv2HucP/93OyK7hOkerhRpkyFfo16EFHz0wlNcTB3CmuJRJs9Zxz7sb2JeVb3dpSjUoG9NPcv+/U+jcpimv3dGfAP8GHU2WavC/SRHhJ70jWPybK3nyum6sTz3B2FdW8Oz8HZw8U2x3eUr5vC0Hc/n52+tpExrE7LsH0Swo0O6SfEqDD/kKQYH+3H9lJ5Y9NpLbB3XgvTVpXPnSMt5aeYCi0jK7y1PKJ20/dIqfzVpHy5BGvH/fYB0q6QYa8hdo3bQxf76hNwt/fQX9o1vy5y92MeblFXy57YheTKWUhXYePs2kWetoFhTI+/cNJqJ5E7tL8kka8tW4rG0zZt89iPfuHkRwI38eTN7ExJlr2KwXUynlsj1H85g0ax1NAv2Ze98QoloGX/pFqk405C/hisvC+eKRy3nx5t6knyjgxtdX89D7mzh4Qi+mUqou9mXlkfjWWgL8hPfvG0J0mAa8O2nI14C/n3DbwGiW/24kj4zqwte7jjHqb9/w/Je7OHW25NI7UEoBcCA7nzv+tQ5wBLwlE3+oi9KQr4WQxgH8dsxlLP/dVYzv156klQcY+dIyZq9Oo6Ss/NI7UKoBS885w53/Wkd5ueH9+wbTuU1Tu0tqEDTk66Bd8yBmTOzL5w+NoHtEKM/M38G1/7eChduP6slZVXvJyRAbC35+jmVyst0VWe7giQLu/Nc6CkvLmHPvYC5r28zukhoMDXkX9IpsTvK9g5k1OQE/P+GBORu5ZeYaUtIuOjfKjzWAP3JVjeRkmDIF0tPBGMdyyhSf+n/gUO5Z7vjXWvIKS5hzz2C6R+gNx+qTJZOGiMjbwDggyxjTy7mtFfAfIBZIA241xlx0aIoVk4bYpbSsnI83ZvLy4u/Jyivi2p5teXxsNzqFX+IracUfeUGlE7nBwT43mbCqRmysI9gvFBPjmG/Uyx3KPcud/1rLifxi5tw7mL4dWthdkk9y+8xQInIFkA+8Vynk/wqcMMa8ICJTgZbGmCcuth9vDvkKBcWlzFqZypsrDnC2pIzbB3bgV6O70KZZNRd5+PgfuboEPz9HC/5CIo6Jpb3YjsOnuOudDZwtLuPduwcRH9PS7pJ8Vr1M/yciscCCSiG/BxhpjDkiIhHAcmNM14vtwxdCvsLx/CJeW7KX5HUZNArw477LOzLlio6ENL5gdhsf/iNXNeCjH/LffJ/Ng3M2EtokkHfuGki3dtpF4052Tf/X1hhzBMC5bFNNcVNEJEVEUrKzs91YTv1q3bQxf5zQi8W/vZKrurbh1SV7ufKl5fx7bfr5I3Giq5l5vrrtyrdMn+7onqssONix3Uv9Z0MGd7+7gQ6tgpn34HANeJvZfuLVGJNkjEkwxiSEh4fbXY7l4lqH8M/EAXz64DA6tg7hD//dzpiXv+GzLYcoLzc++UeuaiEx0XH+JSbG8e0tJsZrz8cYY3h50R6e+GQbwzqF8dEDQ2nXXO9FYzftrqlHxhiW7s7ipa/2sPtoHt3aNeOxa7ty9cavkaemQUaGowU/fbpX/pGrhqu4tJypn37Hp5sOMTE+ir/c1JtAvV1wvblYd407pz+fD0wGXnAuP3PjsbyCiDCqe1uu6tqGz787zMuLv+ee2SnEx8Tx2NKNDOkYZneJStXa6cISHpyziVX7jvOb0ZfxyKjOiIjdZSknSz5qRWQusAboKiKZInIPjnAfIyJ7gTHOdc9Uz+PU/fyECf0i+fq3VzL9xl5knizg9qS1/Ozt9WzLPOXWYytlpSOnznLrzDWsPZDDjIl9+dXoLhrwHsay7hor2NJd4wHj1AtLynhvTRqvL99PbkEJP+ndjt+O6Vrzy76Tk2Gadveo+rXryGnuemcD+UWlvDFpAJd38b1zat6iXoZQWsGWkPegIWynC0t4a2Uqs1Y6xtjf2D+Kh6/uTOzFbuLkAR9SquFZ8X02DyZvomnjAN65a6BexWozDfmL8cBx6jn5Rby+fD9z1qZTWm64sX8kD11VTdh70IeU8n2lZeX8fek+Xlu6l65tm/HOXQN1sg8PoCF/MR4cklmnC5n5zQGS1znC/oZ+kT9u2Xvgh5TyTUdOneVXH2xhfeoJbh4QxZ8m9PzxxX3KFnZdDOUdPHicepvQIJ7+aQ9WPn4VPx8Wy4LvDjPq5W949MOtpB0/43iSXkyl6sHS3cf4yasr2X7oFC/f2pe/3drXuoDXG/S5lYa8F1yM0iY0iD+M68HKJ6oI+6ef99gPKVVP3BiSxaXl/HnBTu5+N4V2zZvw+cMjuGlAlGX7bwh34bSbdtd4oay8QpK+OcCcdekUl5ZzQ/NiHnz/BTpv36CjaxoaN554z8gp4OG5m9iaeYqfDY3h9z/pTlCgv4sFX8CDu0u9ifbJ+6jKYV9YUs6YHm154MpOere/hsRNIbngu8M8+ck2ROCvt/RhbK+IOu/rovSckiW0T97dbOpTbNMsiKfG9eDbJ67mkVFd2JB2gpvfWM2tM9ewbHeWzlLVEGRk1G77JRSWlPHkp9t46P3NdG7blC8eudx9AQ96TqkeaMi7ygP6FMOaNua3Yy7j2yeu5g/jepB5soC73t3Ada+uZN7mTPfPP+vtJ868uX4LQ3Jb5ikm/ONb5q7P4IErO/Hh/UPp0Cr40i90hQcPfPAZxhiP+RcfH2+8TkyMMY54P/9fTIxtJRWXlpmPUw6aMS8vNzFPLDDDnl9i3l51wJwpKrH+YHPmGBMcfP57Dw52bPcGWr/JPVNsnpq3zcROXWDin1tslu0+5saCqzBnjuPvRcSx9JbfvQcBUkw1uap98q7y4D7F8nLDsj1ZzPxmPxvSTtIyOJD/NySGxCExtA216Baw3n7izNvrhzrf1sIYwyebDvH8l7s4WVDMz4bG8ttrLiM0KLBejq+soyde3clLQiIl7QQzvznAkt3H8BdhbK92/HxYLPExLV27oZQHf8jViLfXX0e7j57mD//dzoa0k/SPbsFzE3rRK7J57Xekt9XwCBry7uRl/5On55zh32vS+TDlIKcLS+kREcrPh8Uyvl/7ug2P85IPuWp5e/21lF9UyiuLv+ed1WmEBgUw9bpuTIzvgJ9fHT/oG9jvz1Pp6Bp38oKLqSqLCQvhqXE9WPv7Ufzlxt6UlRse/+Q7hjy/hOf/t4vMkwWX3kll3n7izIr6veDErTHGcRHd35bz1qpUbk2IYumjI7ltYHTdAx4sH92jrKct+QbOGMPaAyd4b00ai3YewxjD6O5tmTwslqEdw2oWAN7eJ+tK/V7wTW5fVj7Pzt/Bqn3H6dk+lOdu6MWAaIuupdCWvEewtbtGRMYCrwL+wFvGmGonD9GQt9fh3LPMWZvOBxsOcuJMMdGtgrl5QBQ3x0cS1dLNQ+m8lQeH3O6jp/nnsv188d1hQhoH8LtrujJpSAz+rrTcL+QFH3INgW0hLyL+wPc4ZobKBDYAdxhjdlb1fA15z1BYUsbC7Uf5aONBVu/PwRgY1imMiQlRjO0ZQZNGFl/a7s088MTt1oO5/GPZPhbvPEZII38mDY3hvss70rppY/cc0Nu/yfkAO0N+KPCsMeZa5/qTAMaY56t6voa858k8WcCnmw7x8cZMMk4U0LRxAOP6RHBLfJTrI3N8gQe15NennuAfy/ax4vtsQoMCuGt4HHcNj6VFcKN6rUPVP7sm8gaIBA5WWs8EBrv5mMpCUS2DeWRUFx6+ujPrU0/w8cZM5m89zAcbDhLXOoRb4qO4oX8kkS0a6MQR06dX3V1RTyeejTGs2nec15buY33qCcJCGvHE2G5MGhJNs5qOd9eWuE9zd0t+InCtMeZe5/r/AwYZYx6u9JwpwBSA6Ojo+PSqWkXKo5wpKuXLbUf4eGMm61JPANA3qjnX9mrHdb0iiLvYdIW+yIaQLCs3LN2dxT+W7WPrwVzahQZx/5UduX1gdO2607RP3Sdod41ym4ycAr7YdoSFO46y9WAuAF3bNmNsr3aM7dWObu2aub9LpwG1RL8/lscnmzL5bPNhjp4upEOrJvziys7cHB9J44AGeJ2DAuwN+QAcJ15HAYdwnHi90xizo6rna8h7t8O5Z1m4/SgLdxxlQ9oJjIHYsOBzLfy+Uc2tD/wG0BLNziti/tbDzNucyfZDp/H3E0ZeFs5NA6K4tmdbAvxduNzFA08cq9qzewjlT4BXcAyhfNsYU21npYa878jOK2LxzmP8b/sR1uzPobTcENE8iBGdWzO8c2uGdQqjjRX3z/HRlmhhSRmLdx7j002ZrNh7nLJyQ+/I5tw0IJKf9m1v3UgZH/39NTR6WwNlq1MFJXy96xiLdx5jzYEcTp0tAaBzm6YM7xTGsM6tGRIXRvPgWt4YC3yqJXq2uIwNaSf44rsjfLntCHlFpUQ0D+KG/pHc1D+SLm2bWX/QBvBNqCHQkFceo6zcsPPwaVbvP863+3PYkHqCsyVl+An0imzOsE6OVn58TMuaTRTtxS3RwpIyNmfksuZADmv2H2fLwVxKygwhjfy5rncENw2IZEhcDa86dkUDOqfhqzTklccqLi1n89sfsnr+ClY3j2VzZDdK/fwRgbjWIfRq35ye7UPp6Vy2DLlgzLcXtUSLS8v5LjOXNftzWL0/h00ZJykqLcdPoHdkc4Z2as3QTmEMim2lF5ypWtGQV57rgpA+ExjE+s4D2Hrvb9jROpYdh05x+FThuae3bx5Ez8gfgr9H+1DaLfgU/6c8pyVqjCE7r4gDx89wIPsMB7Lz2XMsj5S0k5wtKQOgR0QoQzuFMaxTGAPjWtX+Hu5KVaIhrzxXDbpbTpwpZufh0+w4fIrtzmXq8TPnuuID/IR2zYNo36IJkS2a0L6F4+eK9YjmQTW/MKiGSsrKOVNUyuHcQg4czz8X5geOnyE1+wx5RaXnnts4wI+O4U0ZFNuSoZ3CGBwX9uNvJEq5QENeea46njg9U1TK7qOn2XUkj8O5Z53/CjmUe5ajpwspKz9/n82CAmgV0oigAH+CAv1oHOBP40A/ggL9CQr0p3GAH0GBfgQF+BPg78fZ4lLyikrJLywlv6iUvPOWJRSW/Li2yBZN6BgeQsfWIcS1DqFjeFM6hofQvnkT9/erqwbNztsaKHVx0dFVt+QvMRF1SOMA4mNaER/T6kePlZUbsvIKOZxbeO4D4FDuWU6dLaGwpIyi0nIKS8rIKywlO6+IYud6oXNZUlZOcKMAmjYOoFmQ41/rpo2IbR1yblvTxo5/bUIb07F1U+Jah2g/uvJIGvLKXm6494u/nxDRvAkRzZsQH2PRfdOV8lI6M5Syl5fNrKWUt9GWvLJfYqKGulJuoi15pZTyYRrySinlwzTklVLKh2nIK6WUD9OQV0opH6Yhr5RSPkxDXimlfJhLIS8iE0Vkh4iUi0jCBY89KSL7RGSPiFzrWplKKaXqwtWLobYDNwFvVt4oIj2A24GeQHvgaxG5zBhT5uLxlFJK1YJLLXljzC5jzJ4qHpoAfGCMKTLGpAL7gEGuHEsppVTtuatPPhI4WGk907ntR0RkioikiEhKdna2m8pRSqmG6ZLdNSLyNdCuioemGWM+q+5lVWyr8sb1xpgkIAkc95O/VD1KKaVq7pIhb4wZXYf9ZgIdKq1HAYfrsB+llFIucFd3zXzgdhFpLCJxQBdgvZuOpZRSqhquDqG8UUQygaHAFyLyFYAxZgfwIbATWAj8UkfWKKVU/XNpCKUxZh4wr5rHpgN1n95HKaWUy/SKV6WU8mEa8kop5cM05JVSyodpyCtITobYWPDzcyyTk+2uSCllEZ3Iu6FLToYpU6CgwLGenu5YB51cWykfoC35hm7atB8CvkJBgWO7Usrracg3dBkZtduulPIqGvINXXR07bYrpbyKhnxDN306BAefvy042LFdKeX1NOQbusRESEqCmBgQcSyTkvSkq1I+QkfXKEega6gr5ZO0Ja+UUj5MQ14ppXyYhrxSSvkwDXmllPJhGvJKKeXDxBjPmTtbRLKBdBd20Ro4blE5dvKV9wH6XjyRr7wP0PdSIcYYE17VAx4V8q4SkRRjTILddbjKV94H6HvxRL7yPkDfS01od41SSvkwDXmllPJhvhbySXYXYBFfeR+g78UT+cr7AH0vl+RTffJKKaXO52steaWUUpVoyCullA/zqZAXkedE5DsR2SIii0Skvd011ZWIvCQiu53vZ56ItLC7proSkYkiskNEykXE64a7ichYEdkjIvtEZKrd9dSViLwtIlkist3uWlwlIh1EZJmI7HL+v/Uru2uqCxEJEpH1IrLV+T7+aPkxfKlPXkRCjTGnnT8/AvQwxjxgc1l1IiLXAEuNMaUi8iKAMeYJm8uqExHpDpQDbwK/M8ak2FxSjYmIP/A9MAbIBDYAdxhjdtpaWB2IyBVAPvCeMaaX3fW4QkQigAhjzCYRaQZsBG7wtv8uIiJAiDEmX0QCgVXAr4wxa606hk+15CsC3ikE8NpPMGPMImNMqXN1LRBlZz2uMMbsMsbssbuOOhoE7DPGHDDGFAMfABNsrqlOjDErgBN212EFY8wRY8wm5895wC4g0t6qas845DtXA53/LM0tnwp5ABGZLiIHgUTgabvrscjdwP/sLqKBigQOVlrPxAvDxJeJSCzQH1hncyl1IiL+IrIFyAIWG2MsfR9eF/Ii8rWIbK/i3wQAY8w0Y0wHIBl4yN5qL+5S78X5nGlAKY7347Fq8l68lFSxzWu/IfoaEWkKfAL8+oJv8l7DGFNmjOmH49v6IBGxtCvN66b/M8aMruFT3we+AJ5xYzkuudR7EZHJwDhglPHwkye1+O/ibTKBDpXWo4DDNtWiKnH2YX8CJBtjPrW7HlcZY3JFZDkwFrDs5LjXteQvRkS6VFodD+y2qxZXichY4AlgvDGmwO56GrANQBcRiRORRsDtwHyba2rwnCcsZwG7jDEv211PXYlIeMXIORFpAozG4tzytdE1nwBdcYzkSAceMMYcsrequhGRfUBjIMe5aa0XjxS6EXgNCAdygS3GmGttLaoWROQnwCuAP/C2MWa6vRXVjYjMBUbiuKXtMeAZY8wsW4uqIxEZAawEtuH4ewf4vTHmS/uqqj0R6QPMxvH/lh/woTHmT5Yew5dCXiml1Pl8qrtGKaXU+TTklVLKh2nIK6WUD9OQV0opH6Yhr5RSPkxDXimlfJiGvFJK+bD/D3YEE9Xrncp1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xgrid = numpy.linspace(x.min(), x.max(), 30)\n",
    "Xgrid_poly_feat = polynomial_features(xgrid, degree)\n",
    "Xgrid_scaled = min_max_scaler.transform(Xgrid_poly_feat)\n",
    "Xgrid_scaled[:,0] = 1 \n",
    "pyplot.scatter(x, y, c='r', label='true')\n",
    "pyplot.plot(xgrid, Xgrid_scaled@params, label='predicted')\n",
    "pyplot.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Observe underfitting & overfitting\n",
    "\n",
    "In the model above, we just randomly picked a polynomial degree of $3$ for the fitted curve. Is it good enough to model our dataset? Should we try higher-order polynomials?\n",
    "\n",
    "We can repeat our study with different polynomial degrees varying from $1$ to $15$, and see what happens.\n",
    "To faciliate this task, we provide you with a script to train the model and plot these fitted curves interactively using `ipywidget`.\n",
    "\n",
    "Run the cell below and drag the slider to see how the curve changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "71f0d094cc19484483274e0b52147109",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=8, description='degree', max=15, min=1), Output()), _dom_classes=('widge…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append('../scripts/')\n",
    "from plot_helpers import interact_polyreg\n",
    "\n",
    "max_degree = 15\n",
    "interact_polyreg(max_degree, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Underfitting\n",
    "\n",
    "When `degree` is $1$, the straight line clearly fails to capture the underlying relationship between $x$ and $y$. Specifically, a line is too simple to explain how far the data are spread out, namely, the **variance** in the data. We say that the linear model **underfits** the data in this case.\n",
    "\n",
    "Underfitting happens when the model is too naive or the weights need to be trained with more iterations. It is often easy to detect, since a large training error is a good indicator.\n",
    "\n",
    "##### Challenge question\n",
    "\n",
    "> Would having more training data help resolve underfitting?\n",
    "\n",
    "### Overfitting\n",
    "\n",
    "As we increase the polynomial degree, the training error (MAE from the figure title) keeps decreasing. \n",
    "But does it mean that the 15th-order polynomial gives the best fit?\n",
    "Probably not.\n",
    "\n",
    "Drag the slider to the right, and you will find that the curve passes exactly through many points and looks very odd.\n",
    "If we are given new data, this model may not predict well because it fits too closely to the old data.\n",
    "As real-world data tend to be noisy (due to missing and erroneous values), our synthetic data also have noise added.\n",
    "Models with high polynomial degrees are so flexible that they fit the noise rather than the true relationship! \n",
    "In this case, these models are **overfitting**, or have a **high variance**.\n",
    "Overfitting usually happens when the model is overly complicated.\n",
    "The high-order polynomials in our example have too many degrees of freedom (weights) for our data.\n",
    "\n",
    "##### Challenge question\n",
    "\n",
    "> Would having more training data help resolve overfitting?\n",
    "\n",
    "Compared to underfitting, overfitting is in general harder to detect because the training error tends to be small. We will discuss how to identify overfitting later in this module.\n",
    "Now let's focus on how to prevent overfitting using regularization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regularization\n",
    "\n",
    "Regularization is used to reduce or avoid overfitting.\n",
    "The idea is to introduce a term in the loss function that penalizes complicated models. In our polynomial regression model:\n",
    "\n",
    "$$\n",
    "\\hat{y} = w_0 + w_1 x + w_2 x^2 + \\cdots + w_d x^d,\n",
    "$$\n",
    "\n",
    "every polynomial term contributes to the overall complexity of the model.\n",
    "One idea is to add constraints to the magnitudes of these weights.\n",
    "\n",
    "A common approach is to add a regularization term $\\lambda\\sum_{j=1}^d w_j^2$ to the mean-squared error loss:\n",
    "\n",
    "$$\n",
    "L(\\mathbf{w})=\\frac{1}{N} {\\lVert \\mathbf{y} - X\\mathbf{w} \\rVert}^2 + \\lambda \\sum_{j=1}^d w_j^2\n",
    "$$\n",
    "\n",
    "This new loss function favors smaller weights, because larger weights will increase the term on the right.\n",
    "As weights are smaller, the model is less likely to overfit by large amplitude higher-order polynomial terms.\n",
    "\n",
    "Notice that we don't penalize the intercept $w_0$ in the regularization term. This can be justified mathematically, but one way to think about it is that a constant term won't affect the overall model complexity.\n",
    "\n",
    "Above, $\\lambda$ denotes the regularization parameter, or the strength of penalty. It controls the tradeoff between fitting the data well (minimizing the first term) and keeping the model simple to avoid overfitting (minimizing the second term). When $\\lambda\\rightarrow 0$, the loss function falls back to standard mean-square error loss. And when $\\lambda$ is large, all penalized weights will be close to 0 after training and our model would approach the constant value $w_0$.\n",
    "\n",
    "##### Note\n",
    "\n",
    "> This regularization is also called $L^2$-penalty, or Tikhonov regularization, while the name _ridge regression_ is also used for this model-building method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's code the regularized mean-squared error loss and set $\\lambda = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regularized_loss(params, model, X, y, _lambda=1.0):\n",
    "    '''\n",
    "    The mean squared error loss function with an L2 penalty.\n",
    "    Arguments:\n",
    "      params: 1D array of weights for the linear model\n",
    "      model : function for the linear regression model\n",
    "      X     : 2D array of input values\n",
    "      y     : 1D array of predicted values\n",
    "      _lambda: regularization parameter, default 1.0\n",
    "    Returns:\n",
    "      float, regularized mean squared error\n",
    "    '''\n",
    "    y_pred = model(params, X)\n",
    "    return numpy.mean( numpy.sum((y-y_pred)**2) ) + _lambda * numpy.sum( params[1:]**2 )\n",
    "\n",
    "gradient = grad(regularized_loss)                                                        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_regularization_params = params.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And train the 3rd-degree polynomial model using gradient descent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration    0, loss = 5434.768, mae = 11.05718775676392\n",
      "iteration  100, loss = 1785.985, mae = 6.983886996350374\n",
      "iteration  200, loss = 1764.260, mae = 6.966518785947395\n",
      "iteration  300, loss = 1763.570, mae = 6.965728405272657\n"
     ]
    }
   ],
   "source": [
    "max_iter = 3000\n",
    "alpha = 0.01\n",
    "params = numpy.zeros(X_scaled.shape[1])\n",
    "descent = numpy.ones(X_scaled.shape[1])\n",
    "i = 0\n",
    "\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "while numpy.linalg.norm(descent) > 0.01 and i < max_iter:\n",
    "    descent = gradient(params, linear_regression, X_scaled, y)\n",
    "    params = params - descent * alpha\n",
    "    loss = mse_loss(params, linear_regression, X_scaled, y)\n",
    "    mae = mean_absolute_error(y, X_scaled@params)\n",
    "    if i%100 == 0:\n",
    "        print(f\"iteration {i:4}, {loss = :.3f}, {mae = }\")\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the optimal weights before and after regularization. We can reuse the `xgrid` to plot both curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "weights without regularization\n",
      "[-22.51572398   6.75930601  41.30788709  30.0105898 ]\n",
      "weights with regularization\n",
      "[-11.13750882  12.48522096  28.26626633  11.09211867]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6bUlEQVR4nO3deVhV1frA8e9iEAREFFBxAJxREec5p7S0MjNNLanMW5nZdKtbWnbLunnr+rPbZDczs6zIMss0K3PKecQpZ3MARBARRUGmczjr98cGRGXm4Bl8P89zns3e5+y93yP4nnXWXvtdSmuNEEII5+Ri6wCEEEJUHUnyQgjhxCTJCyGEE5MkL4QQTkySvBBCODE3WwdQWEBAgA4NDbV1GEII4VB27NhxVmsdWNRzdpXkQ0NDiY6OtnUYQgjhUJRSscU9J901QgjhxCTJCyGEE5MkL4QQTsyu+uSLYjKZiI+PJysry9ahCBvx9PSkYcOGuLu72zoUIRyO3Sf5+Ph4atSoQWhoKEopW4cjrjOtNSkpKcTHx9O4cWNbhyOEw7H77pqsrCz8/f0lwd+glFL4+/vLNznhvKKiIDQUXFyMZVSUVQ9v9y15QBL8DU5+/8JpRUXB+PGQkWGsx8Ya6wCRkVY5hd235IUQwmlNmXI5wefLyDC2W4kkeSt56623iLLy16yKeOihh1i4cGG59pk1axZffvlluc+1Zs0aNm3aVOnjCHHDiosr3/YKcIjuGkewfPlyFixYUO79zGYzbm62+zWYzWYmTJhQoX3XrFmDj48PPXv2BKjwcYS4YQUHG100RW23EmnJl2L69Ol88MEHADz77LPcfPPNAKxatYr7778fgIsXL5KTk0NgYCCxsbEMGDCAiIgIBgwYQFwRn8hTp05l/Pjx3HrrrTz44IMkJyczYsQIunTpQpcuXdi4cSMAycnJ3HLLLXTs2JHHHnuMkJAQzp49S0xMDOHh4QXHmzFjBlOnTr3mPG+88QZdunQhPDyc8ePHkz8LWL9+/Xj55Zfp27cv77//PlOnTmXGjBkkJCTQvn37goerqyuxsbH8/PPPdOvWjQ4dOjBw4ECSkpKIiYlh1qxZvPvuu7Rv357169cXHAdg9+7ddO/enYiICO6++27Onz9fcO5JkybRtWtXWrRowfr16630mxLCAU2bBl5eV27z8jK2W4lDteRf/3k/BxIuWvWYrev78tqdbYp9vk+fPrzzzjs8/fTTREdHk52djclkYsOGDfTu3RuAlStXMmDAAACefPJJHnzwQcaOHcvcuXN5+umn+emnn6457o4dO9iwYQPVq1dnzJgxPPvss9x0003ExcUxaNAgDh48yOuvv87NN9/MSy+9xLJly5g9e3a53tuTTz7Jq6++CsADDzzA0qVLufPOOwFITU1l7dq1AAUfEPXr12f37t0AfPTRR6xdu5aQkBB8fX3ZsmULSinmzJnD9OnTeeedd5gwYQI+Pj784x//AIwPvnwPPvggH374IX379uXVV1/l9ddf57333gOMbw/btm3j119/5fXXX2flypXlel9COI38i6tTphhdNMHBRoK30kVXcLAkbwudOnVix44dpKWl4eHhQceOHYmOjmb9+vUFLfxly5Yxbtw4ADZv3syPP/4IGIn1xRdfLPK4Q4cOpXr16oDxIXHgwIGC5y5evEhaWhobNmxg0aJFAAwePJhatWqVK/Y//viD6dOnk5GRwblz52jTpk1Bkh89enSx+23cuJE5c+YUtLLj4+MZPXo0iYmJ5OTklDpe/cKFC6SmptK3b18Axo4dy8iRIwueHz58OGD828bExJTrPQnhdCIjrZrUr+ZQSb6kFndVcXd3JzQ0lM8//5yePXsSERHBH3/8wbFjx2jVqhUA27Zt4+OPPy5y/+KG/3l7exf8bLFY2Lx5c0HSz1fcJOtubm5YLJaC9aLGkGdlZTFx4kSio6Np1KgRU6dOveJ1hc9fWGJiIg8//DBLlizBx8cHgKeeeornnnuOoUOHsmbNmiK7hsrDw8MDAFdXV8xmc6WOJYQomfTJl0GfPn2YMWMGffr0oXfv3syaNYv27dujlGL//v2EhYXh6uoKQM+ePfn2228BiIqK4qabbir1+LfeeiszZ84sWM/vMrnpppsKLuYuX768oF+7bt26nDlzhpSUFLKzs1m6dOk1x8xP6AEBAaSnp5dpxI3JZGLUqFH85z//oUWLFgXbL1y4QIMGDQCYN29ewfYaNWqQlpZ2zXFq1qxJrVq1Cr4JfPXVVwWteiHE9SVJvgx69+5NYmIiPXr0oG7dunh6ehb0x//2228MHjy44LUffPABn3/+OREREXz11Ve8//77pR7/gw8+IDo6moiICFq3bs2sWbMAeO2111i+fDkdO3bkt99+IygoiBo1auDu7s6rr75Kt27dGDJkCGFhYdcc08/Pj0cffZS2bdsybNgwunTpUmocmzZtYvv27bz22msFF18TEhKYOnUqI0eOpHfv3gQEBBS8/s4772TRokUFF14LmzdvHi+88AIRERHs3r274NqAEOL6UsV1CdhC586d9dWThhw8eLCgW8Qe3XLLLXz55ZcEBQVZ/djZ2dm4urri5ubG5s2befzxxwta+Tcae/87EMKWlFI7tNadi3rOofrk7dGKFSuq7NhxcXGMGjUKi8VCtWrV+PTTT6vsXEII52SVJK+U8gPmAOGABv4GHAa+A0KBGGCU1vq8Nc53o2jevDm7du2ydRhCCAdmrT7594FlWuswoB1wEJgMrNJaNwdW5a0LIYS4jiqd5JVSvkAf4DMArXWO1joVuAvIH4oxDxhW2XMJIYQoH2u05JsAycDnSqldSqk5SilvoK7WOhEgb1mnqJ2VUuOVUtFKqejk5GQrhCOEECKfNZK8G9AR+Fhr3QG4RDm6ZrTWs7XWnbXWnQMDA60QjhBCiHzWSPLxQLzWemve+kKMpJ+klAoCyFuescK57NLtt99Oamoqqamp/O9//yvYvmbNGoYMGWLDyCA0NJSzZ8+Wa59HHnnkijILZfXFF1+QkJBQ6eMIIayn0klea30aOKmUapm3aQBwAFgCjM3bNhZYXNlz2atff/0VPz+/a5J8Rdj6Nv/c3FzmzJlD69aty73v1Um+oscRQliPtUbXPAVEKaX+BNoD/wbeBm5RSv0F3JK37nDKUmo4v7U8efJkjh07Rvv27XnhhRcASE9P55577iEsLIzIyMgi69FcXfp3x44d9O3bl06dOjFo0CASExMB2L59OxEREfTo0YMXXnihoNzwF198wZNPPllwvCFDhrBmzZprzjNs2DA6depEmzZtrqho6ePjU3AH7ebNm+nXrx/R0dEsWbKk4M7Xli1bFhQmK6qE8cKFC4mOjiYyMpL27duTmZlZcByA+fPn07ZtW8LDw5k0adIV554yZQrt2rWje/fuJCUlVewXJYQoklXGyWutdwNF3W01wBrHL/DbZDi916qHpF5buK34z5+ylBrO9/bbb7Nv376Cu1LXrFnDrl272L9/P/Xr16dXr15s3LixyHo2+aV/TSYTffv2ZfHixQQGBvLdd98xZcoU5s6dy7hx45g9ezY9e/Zk8uTyj0idO3cutWvXJjMzky5dujBixAj8/f25dOkS4eHhvPHGG1e8fujQoQwdOhSAUaNGFdSfKaqE8T333MPMmTOZMWMGnTtf+aeQkJDApEmT2LFjB7Vq1eLWW2/lp59+YtiwYVy6dInu3bszbdo0XnzxRT799FNeeeWVcr83IUTRpHZNKa4uNdyjR4+CUsNXJ/midO3alYYNG+Li4kL79u2LLa2bX/r38OHD7Nu3j1tuuYX27dvz5ptvEh8fT2pqKmlpaQWzMI0ZM6bc7+WDDz4oaDGfPHmSv/76CzCqQY4YMaLY/aZPn0716tV54oknAKOEcbdu3Wjbti2rV69m//79JZ53+/bt9OvXj8DAQNzc3IiMjGTdunUAVKtWreC6hZQeFsL6HKusQQkt7qpSllLDJckvqwsll9bNL/2rtaZNmzZs3rz5iufzK1AWpSylh9esWcPKlSvZvHkzXl5e9OvXr+B1np6eBVU0r7Zq1Sq+//77gqRcWgnjopRUH8nd3b2gHLOUHhbC+qQlXwYllRourLjSu+XRsmVLkpOTC5K8yWRi//791KpVixo1arBlyxaAgnLGYFwT2L17NxaLhZMnT7Jt27ZrjnvhwgVq1aqFl5cXhw4dKjhOSWJjY5k4cSILFiwoqHVfUgnj4t5/t27dWLt2LWfPniU3N5f58+dL6WEhrhPHasnbSO/evZk2bRo9evTA29v7ilLDhfn7+9OrVy/Cw8O57bbbuOOOO8p9rmrVqrFw4UKefvppLly4gNls5u9//ztt2rThs88+49FHH8Xb25t+/fpRs2ZNAHr16kXjxo0LLmx27NjxmuMOHjyYWbNmERERQcuWLenevXupsXzxxRekpKRw9913A8b0gL/++mtBCePQ0NArShg/9NBDTJgwgerVq1/xTSQoKIi33nqL/v37o7Xm9ttv56677ir3v40Qovyk1LADSU9PL5it6e233yYxMbFM9eqdgfwdCGf2/II9NAn05on+zSq0f0mlhqW7xoH88ssvtG/fnvDwcNavXy+jUIRwAodOX+SHnfFkmy2lv7gCpLvGgYwePbrECbiFEI7nw9VH8a7myt96hVbJ8aUlL4QQNnL0TBq/7k1kbM9Q/LyqVck5JMkLIYSNzFx9FE83Vx6+qXGVnUOSvBBC2MCJs5dYsieBB3qE4O/jUfoOFSRJXgghbOCjP47i7urCo72bVOl5JMmXgTWqSwohRL6T5zJYtOsUY7oFE1ij6lrxIEm+TIpL8rm5uTaIRgjh6P635iiuLooJfZtW+bmcL8lHRUFoKLi4GMuoqEofsnAJ4S5dutC/f3/GjBlD27ZtiYmJKSj5CzBjxgymTp0KwLFjxxg8eDCdOnWid+/eHDp0qNKxCCEc26nUTBbuiGd050bU9fWs8vM51zj5qCgYPx4yMoz12FhjHSAyssKHLVxCeM2aNdxxxx3s27ePxo0bl1g1cfz48cyaNYvmzZuzdetWJk6cyOrVqyschxDC8c1acwyACf2qvhUPzpbkp0y5nODzZWQY2yuR5K/WtWvXggk0ipOens6mTZsYOXJkwbbs7GyrxSCEcDynL2Tx3faT3NOpIQ38ql+Xc1olySulYoA0IBcwa607K6VqA98BoUAMMEprXXy9XGuIiyvf9grKLwsMxZf5tVgs+Pn5FUwgIoQQn6w7Rq7WTOxXsRo1FWHNPvn+Wuv2hYrkTAZWaa2bA6vy1qtWcHD5tpdRSSWE69aty5kzZ0hJSSE7O5ulS5cC4OvrS+PGjfn+++8Bo6b6nj17KhWHEMJxnUnL4putcdzdoQGNantdt/NW5YXXu4B5eT/PA4ZV4bkM06aB11X/eF5exvZKKFxCOH/u1nzu7u4F86MOGTKEsLCwgueioqL47LPPaNeuHW3atGHxYqedy1wIUYo5609gyrVUuNJkRVml1LBS6gRwHtDAJ1rr2UqpVK21X6HXnNda1ypi3/HAeIDg4OBOsbGxVzxf7hKzUVFGH3xcnNGCnzbNqv3xwjak1LBwZCnp2dz0nz8YHF6Pd0e3t/rxSyo1bK0Lr7201glKqTrACqVUmccKaq1nA7PBqCdf6UgiIyWpCyHsymcbTpBlzr3urXiwUneN1johb3kGWAR0BZKUUkEAecsz1jiXEEI4ktSMHL7cHMsdbYNoVsfnup+/0kleKeWtlKqR/zNwK7APWAKMzXvZWKDCHdL2NHuVuP7k9y8c2dyNMaRnm3ny5uvfigfrdNfUBRblTWrtBnyjtV6mlNoOLFBKPQzEASNLOEaxPD09SUlJwd/f/5qJs4Xz01qTkpKCp2fV3xkohLVdzDLx+cYTDG5Tj7B6vjaJodJJXmt9HGhXxPYUYEBlj9+wYUPi4+NJTk6u7KGEg/L09KRhw4a2DkOIcpu3MYa0LNu14sEB7nh1d3cv9e5SIYSwN+nZZj7beIKBreoQ3qCmzeJwvgJlQghhBz7fcILUDBNP3dzcpnFIkhdCCCtLSc/mk3XHGdSmLu0a+dk0FknyQghRGUWUN5/5x1EyTbm8ODistL2rnCR5IYSoqPzy5rGxoDXExhL3j1f4etMJRnVuRNPA6z8u/mqS5IUQoqKKKG/+Tud7cDWZ+PtA2/bF55MkL4QQFXVVGfN9dZqwuE0/Ht7+03WZ9aksJMkLIURFXVXG/D/9HqJWxgUeS9hqo4CuJUleCCEqqlB58w0h7VjfuCNPRi/Cd+o/bRzYZXZ/M5QQQtitvIq3limv8Ha/cTS4lML9jw+zq0q4kuSFEKIyIiNZGt6fffN38e7odnh0sK8SHNJdI4QQlZBjtjDj98O0CvLlrnYNbB3ONSTJCyFEJczfFkfcuQwm3xaGi4v9VcqVJC+EEBWUlmXig1V/0bOpP32aB9g6nCJJkhdCiAr6dP0JUi7lMPm2MLud78JqSV4p5aqU2qWUWpq3XlsptUIp9Vfe8ppJvIUQwlGdSctizvrjDIkIIqKhX+UOlnkeLBarxHU1a7bknwEOFlqfDKzSWjcHVuWtCyGEU/hw1VFyzBb+cWvLyh3IlAVfDYefHrdOYFexSpJXSjUE7gDmFNp8FzAv7+d5wDBrnEsIIWztxNlLzN8Wx5huwYQGeFf8QFrD0mchYSe0Hmq9AAuxVkv+PeBFoPD3jbpa60SAvGUdK51LCCFsasbvh/Fwc6n8hCDbZsOeb6DfSxB2h3WCu0qlk7xSaghwRmu9o4L7j1dKRSulomUeVyGEvdt9MpVf9ibyaJ8mBNbwqPiBTqyHZS9Byzugz4vWC/Aq1mjJ9wKGKqVigG+Bm5VSXwNJSqkggLzlmaJ21lrP1lp31lp3DgwMtEI4QghRNbTWvPXrQQJ8qvFI7yYVP1BqHHw/Fvybwt2zjAlHqkilj6y1fklr3VBrHQrcC6zWWt8PLAHG5r1sLLC4sucSQghb+vnPRLaeOMezt7TAx6OCVWFyMuDbSMg1wb3zwdPXukFepSpr17wNLFBKPQzEASOr8FxCCFGl0rPNTPvlAG0b1OTeLsGl71AUreHnZ+D0XhjzHQQ0s26QRbBqktdarwHW5P2cAgyw5vGFEMJWPlz1F0kXs5l1fydcK1q+YPNHsHcB3PwKtBhk3QCLIXe8CiFEKY6eSeOzDScY3bkRHYIreF/nsT9gxT+h1VDo/Q/rBlgCSfJCCFECrTWvLt6PVzVXXhxcwRufzp2AheMgMAyGfQzXsQSCJHkhhCjBL3sT2XQshRcGtcTfpwJDJnMuGRdatQXujQIPH+sHWQKZNEQIIYpxKdvMm0sP0qa+L2O6hZT/AFrD4icg+SBEfg+1KzHssoKkJS+EEMX4cPVRTl/M4o27wit2sXXje7B/EQx4DZoNtHp8ZSFJXgghinAsOZ3PNhxnZKeGdAqpwMXWw8tg5evQZjj0esb6AZaRJHkhhLiK1pqpS/ZT3d2VSbeFlf8ACbuMC61BEXDXzOt6ofVqkuSFEOIqv+07zfq/zvL8rS0JKO/F1vOx8M1o8PKHMQugWiWqVFqBXHgVQohCMnLMvLn0AK2CfInsVs47WzPPQ9RIo0b8g0ugRr2qCbIcJMkLIUQhM1cfJeFCFh/c1wE313J0dpiz4bsH4NxxeOBHqFOBbp4qIEleCCHyHE9O59P1xxnesQGdQ2uXfUetYfGTELMehn8KjftUXZDl5Bx98lFREBpqlOsMDTXWhRCiHLTWvLZkP55urrx0W6vy7bz6zcs1aSJGVU2AFeT4LfmoKBg/HjIyjPXYWGMdIDLSdnEJIRzK7/uNi62vDmldvslAdsyD9TOg44PXtSZNWTl+S37KlMsJPl9GhrFdCCHKIDMnl38tPUhYvRo82KMcd7YeXWnM0dp0ANzxX5sOlSyO47fk4+LKt10IIa7yzvLDnErN5Lvx3ct+sTXxT1gwFuq2hlHzwNW9aoOsIMdvyQcXM8SpuO1CCFHIjthzfLbxBJHdgunWxL9sO12Ih29GgWdNGPM9eNSo2iArwRoTeXsqpbYppfYopfYrpV7P215bKbVCKfVX3rKCRZhLMW0aeHlduc3Ly9guhBAlyDLl8sL3f1K/ZnVeur2MF1uzLkDUKKO6ZOT34BtUtUFWkjVa8tnAzVrrdkB7YLBSqjswGViltW4OrMpbt77ISJg9G0JCjP6wkBBjXS66CiFKMeP3wxw/e4np90SUbc5WU5YxFv7sYRj9FdRtU/VBVlKl++S11hpIz1t1z3to4C6gX972eRjTAk6q7PmKFBkpSV0IUS7RMZe7aXo1Cyh9h1wTfP8QnFgLd38CTfpVdYhWYZU+eaWUq1JqN3AGWKG13grU1VonAuQt6xSz73ilVLRSKjo5Odka4QghRImyTLm8sLAc3TSWXPhxPBz5De54B9rdW/VBWolVkrzWOldr3R5oCHRVSoWXY9/ZWuvOWuvOgYGB1ghHCCFKNOP3w5w4e4n/K0s3jcUCPz8N+3+EW96ALo9cnyCtxKqja7TWqRjdMoOBJKVUEEDe8ow1zyWEEBVRuJumZ2ndNFrD7y/Brq+h7ySb1oWvKGuMrglUSvnl/VwdGAgcApYAY/NeNhZYXNlzCSFEZWTmlLObZvW/YOss6PEk9Hup6gOsAta4GSoImKeUcsX40FigtV6qlNoMLFBKPQzEASOtcC4hhKiwGcuNbppvHulWejfNuhmw/h3o9BDc+qZd3s1aFtYYXfMn0KGI7SnAgMoeXwghrCE65hxzN57g/u5l6KbZMstoxbcdZbflCsrK8e94FUKIUhTupplcWoXJnV/BskkQNgSGfQwurtcnyCri+LVrhBCiFGXuptm7EJY8ZRQcu2cuuDp+ipSWvBDCqW0vazfNoV9h0WMQ3ANGfw1u5Zzb1U5JkhdCOK3MnFxeXPgnDfyqlzwRyF8r4PuxUC8CxnwH1byKf62DkSQvhHBa/5d309P0ERF4F9dNc2AJzL8PAsPg/h/A0/f6BlnFJMkLIZzSqoNJzN14ggd7hBTfTfPnAqMeTf0OMPZn8CrHvK4OQpK8EMLpnErN5Pnv99AqyJeXi7vpaccXRj2akJ7wwCKo7nc9Q7xuJMkLIZyKKdfCU9/sxGS28L/Ijni6FzEEcsvH8PMz0GygURPew+f6B3qdOP74ICGEKGTG74fZGZfKB/d1oHGA97UvWDfDuNGp1Z0w4jOnGUVTHEnyQginsepgEp+sO05kt2CGtqt/5ZNaG8l9/TvGnazDPnaKcfClcf53KIS4IeT3w7cO8uWfQ1pf+aTWsOwl2PoxdBwLQ951+DtZy0qSvBDC4eX3w5tzNR9d3Q9vyYWlz8LOedDtcRj8lkPXoikvufAqhHBsUVHMGP4cO+NSeWv1JzT+/afLz+WaYdEEI8H3fv6GS/AgLXkhhCOLimLVv2fxyZ2TuX/nL9y5YRHs/N14btRw+OEROLQUbv4n9PmHbWO1EWnJCyEc1qlpM3h+4ERaJx3jldVzjI0ZGfCvl+CLIXDoF7hteskJPioKQkPBxcVYRkVdj9CvG2nJCyEckinXwlOdIjG7uPHR4v/gmWsynvB3gUHnIWk/jP7KGCpZnKgoGD/e+GAAiI011gEiI6v2DVwn1pj+r5FS6g+l1EGl1H6l1DN522srpVYopf7KW9aqfLhCCGGY8fthdjZoxVvLPqTx+QRjY7ArPOwNnq7w0C8lJ3iAKVMuJ/h8GRnGdidhje4aM/C81roV0B14QinVGpgMrNJaNwdW5a0LIUSl5Y+Hv792NnfG7TA2tnWDB7wgA2j1OjTsVPqB4uLKt90BVTrJa60TtdY7835OAw4CDYC7gHl5L5sHDKvsuYQQIv58RsF4+FeevQs++QTurAPDvSDZHbq8Aw89U7aDBQeXb7sDsuqFV6VUKMZ8r1uBulrrRDA+CIA6xewzXikVrZSKTk5OtmY4Qggnk5Zl4uEvosm15I2Hd7GAz2bomAUR98JHifDgo2U/4LRp4HVV7XgvL2O7k7BakldK+QA/AH/XWl8s635a69la685a686BgYHWCkcI4WTMuRae/GYXx5LTmXV/Jxp7m+DrEbD7a+g7Ge6eBW7VynfQyEiYPRtCQozx8yEhxrqTXHQFKyV5pZQ7RoKP0lr/mLc5SSkVlPd8EHDGGucqSkp6NiM+3sSaw1V2CiGEDWmtef3nA6w9ksy/hoXTKyAD5g6G2I1GDZr+L1X8JqfISIiJAYvFWDpRggfrjK5RwGfAQa31fws9tQQYm/fzWGBxZc9VnMQLWaSkZ/PQ59t5cO42Dp9Oq6pTCSFs4ItNMXy1JZbH+jThvqAkmDMQLibA/T9C+zG2Ds+uKa115Q6g1E3AemAvYMnb/DJGv/wCIBiIA0Zqrc+VdKzOnTvr6OjoCsWRY7bw1ZZY3l95hPRsM/d2DebZgS0IrOHcZUSFcHarDibx6JfRDAyrw6ywXbj8/jL41jfmYq1TwrytNxCl1A6tdecin6tskremyiT5fOcv5fD+qr/4ekssnu6uTOzflL/1alz0xAFCCLt2IOEi98zaROsAN76t/y1u+76H5rfC8NlQXW69yVdSkneOsgbZ6bDwYUj8k1re1Zg6tA2/P9uH7k38mb7sMAPeWcvPexKwpw80IUTJki5m8fC87bTxPMt8l1dw27cQ+k+B+76TBF8OzpHkzxyAoyvhk95GQaJzJ2ga6MOcsZ355pFu+FZ356n5uxj+8SZ2xp23dbRCiFJk5Jh5ZF40nTI38y0v4X7pNNy/EPq+aNSYEWXmPN01mamw8X1j7kaLGTqPgz4vgk8guRbNDzvi+b/lh0lOy+bOdvWZNLglDWt5lXpYIcT1ZbFoHv9qG+3+mslEtyUQ1B5GfQm1Qmwdmt26ofrkuZgIa9+GnV+Bmyf0fAp6PgkeNbiUbWbW2mN8uv44Fg2P3NSYif2b4eMhddqEsBfvLd5E5+h/cJPrfmMWp9umg7unrcOyazdWks939i9jPscDi8ErAPq8YLTu3TxISM3k/34/zKJdpwjwqcbzt7ZkVOdGuLrcWJMJCGFvli9fSvjGpwh0Scd96LvQ4X5bh+QQbswkny9+B6x8DWLWg1+wMXlA+D3g4sLuk6m8ufQA0bHnCatXg1fuaM1NzQOse34hROm05tgv79Fo+7+44BZIrXHf4tawg62jchg3dpIHYxLfY6tgxVRI2gt12hh3yIUNQQO/7TvNW78d5OS5TAaE1eGl21vRrI6P9eMQQlwr7TQXvn2MmqfWsM2tM62emE+NWkWWuhLFkCSfz2KBfT/Amrfg3DGoFwH9X4YWg8kyW5i3KYaZq4+Sacrl/u4hPDOgObW8y1kLQwhRdvsXYV7yd0xZGXziOY77Jr5O3ZrVbR2Vw5Ekf7VcM+xdAGv/A+djoH4HY/xts4GcvZTDeyuP8M3WOHw83Hh6QHMe6BGCh5vcTCWE1WSmwm8vwp/fsZemvO35HDMeH0GQJPgKkSRfnFwT7JkPa/8PLsRBwy5Gy75Jf46cSWfaLwdZeySZ4NpeTL4tjNvC66FusJnehbC642vhp4notERmMYL51UbyzYSbZEhzJUiSL405xyhXum4GXDwFwT2MZN+4D+uOJPPvXw9y6HQanUNqMeWOVnQIlrvthCg3UyasegO2/I9sv6Y8fPFRjrm34LvxPQj2lwRfGZLky8qcDTu/hPXvQFoihPaGPi+QG9KbhTvjmbH8CMlp2QyJCGLS4DAa1ZY/TCHKJGE3/Dgezh4mte1D3HHgFkwunnz3WA8aB3jbOjqHJ0m+vExZsOML2PBfSE+CBp2h9/NcCh3IJ+tjmL3uGBYLjOsVysT+zahZ3d3WEQthn3LNsPFdWPM2eAeS2P+/DP2tGlrDt+O7yyg2K5EkX1GmLNgdBRvfg9Q4qNMabnqW041uZ8bKY/ywMx6/6u78fWALxnQLxt1VamoIUSBuC/zyPCTtg/ARnOzxBiPnHcaUa2H++O60qFvD1hE6DeevQllV3D2hy8Pw1C64ezZoC/z4KPW+7MWMxjtYOqEzrYJ8eW3Jfga9u45l+05LpUtRflFREBpqFN4KDTXWHdmlFFj8BMwdZIyiGf01J2+eyb1fHSHLnMvXj3STBH8dSZIvC1c3aDcaHt8Mo6Ogem1Y+ixtvu9DVOttfDGmFS4uiglf7+CeWZuJjilxbpRrOdt/clF2UVEwfjzExho37cXGGuuO+DdgscCOeTCzE+z5Fno9A09s5VTQQO77dAtpWSa+frgbrYJ8bR3pDcUq3TVKqbnAEOCM1jo8b1tt4DsgFIgBRmmtS6zza3fdNcXRGk6sNS7QnlgH1Wth6fwoi6vdxlvrznEmLZtBbery4uAwmgaW0ueY/588I+PyNi8vp5tMWBQjNNRI7FcLCTHmG3UUiX/CL89B/HYI6QV3vAN1WnEqNZMxn27hXHoOXz/SjXaN/GwdqVOq8j55pVQfIB34slCSnw6c01q/rZSaDNTSWk8q6TgOk+QLO7nduEB7+FdwrYa5zT0scBvKv3e4kGnK5d4ujXhmYHPq1Cimip6z/CcXFePiYjQarqaU0TK2d1kX4Y9/w7ZPjG+4t74J7e4FpdifcIFxn28nMyeXL/7WlU4hMvS4qlyXC69KqVBgaaEkfxjop7VOVEoFAWu01i1LOoZDJvl8Z4/C1o9h9zdgyiAnuDcL3Iby+qH6uLu58WjvJozv0wTvq8saO/p/clE5jvohrzXs/xGWvWyMQOv8Nxjwz4IZm9YeSWbi1zvwre7O5+O6EFZPumiqkq2SfKrW2q/Q8+e11td8lCulxgPjAYKDgzvFFvUH70gyzsHOebB1NqQlkOPXlB/ch/L6yQh8fHx5ZmBz7u3S6PJIHEf9Ty6swxG76wpXdg1qD0P+Cw06FTz93fY4Xl60j+Z1fPhiXFfq1ZRa8FXNrpN8YQ7dkr9argn2/wRbPoKEXZg9/PjZbRBvpfTBy78Bz97Sgjsj6uMy/xvH+08urCsqCqZMgbg4CA6GadPs83d/5pAxR8OhpcYcDf0mGy14F6Ouk9aad1cc4YPVR+ndPID/RXakhqfcQ3I9SHeNLWltjBfe8hH60C9oXFjv1p1Zl/pyPrAbLwwO4+YdK1GvOMB/cnFjSo0zbmbaMx/cvY3Z1npMBI/LwyBzzBYm//gnP+48xchODfn38LZy38h1VFKSr8p575YAY4G385aLq/Bc9kspCOkBIT1Q506gts+hz+4o+lbbQNzF+sz7uj9f1h/K46t30L2Jv62jFeKy9GRjBFn0Z4CC7hPhpufA+8q/04tZJiZ+vZMNR8/y7MAWPD2gmRTysyPWGl0zH+gHBABJwGvAT8ACIBiIA0ZqrUscQG6zlvz1/rpsyoIDi7FEf4bLya1k487S3G4cqD+CYUOG01aGmQlbyroIm2fC5o/AlAHtI42umZoNr3lp4oVMxn2+naNn0nl7RAT3dLr2NaLqSVmDktj6wlfSfszb5mLZ/S3VctM5ZGnE7rp303Xo4zRpVL9sx3CUPl1h33IyIHqu0XrPPAet74L+r0BgiyJffjDxIuM+3056tpmP7+9I7+aB1zlgkU+SfEnsZXRLziUydy3gwvpPqJd+kAztwb5aAwju9zfqRQwwhloWxdYfUsLxpZ+BbbNh+2dGcm/SHwa8Cg06FrvLuiPJTIzaiY+HG5+P6yJ3sdqYJPmS2OE49dSjWzm+bCYtklfgozI5714X1e5e/Lo/AAHNr3yxvXxICceTfNjoltnzHeTmQMvbjYuqIT2K3cWca+GD1Uf5cPVftKxbg8/HdZHZnOyAJPmS2HGSPJNyjvVLvyTw2I/0Un/iqjTZdTvg0SkS2gw3LoDZ4YeUsGNaQ8wG2PQh/PU7uHlC+zHQ/QkIaFbirokXMnnm291sO3GOER0b8sZdba69uU/YhCT5kjhAd8eZi1lErdxG1q4F3KXW0dolFu3ijmoxCP73B2xKhNyrdrKDDylhR3JNcGAxbPoAEvcY49y7jjeqrHoHlLr76kNJPL9gD9lmC28OC2d4RyteYJVrSpUmSb40DvJHdiYti0/WHmf7lnUMVesY5bEZX/M5yNRwyAQHzXDcDB729SElqlhJf79pScb49u1z4MJJ8G8OPZ4w6su4l97NkmO2MH3ZIeZsOEGrIF9mjulQetG98sZu540sRyBJ3smcScti9trjzN96nM65f/KU5zo6ZG/F1d0COYB/F7jtaWg2EKrJFIVOragkWcMLpk+AGifh6ErQuUZlyJ5PQfNBxV/Ev0pcSgZPzd/JnvgLPNgjhJdvb4Wnu6t147fj7lJHIkneSeUn+6+3xmIxZfNE6Ckia+wmIH6lMUrC3Qua3wKthkKLQVfcoSicROEkWc8FOlSDcDfwcoEa9Y0We/sx116wL8XSPxN46Ye9KAXT74lgcHiQ9WMHuaZkJZLkq5qNu3tS0rOZtzmWLzfHkJphontITSa3Oku79HWoQ0uNKoGuHtD0Zgi7A5oNAN8yjsEX9s3HFdq6QXt3qOsKZg2HzLDHBIezCurKlFWWKZfXfz7A/G1xdAj244N7O1TthPXSkrcKmf6vKtnBzD7+Ph48d0sLNk66mX8OaU1cajbDlrlz29FhLLp5BeaxvxoX2E7vhSVPwn9bwf96wvJ/wvG1YM6uXACOPrOVo8V/KQV2z4f598FzPjDIE0zA0kx4Jw1+yARTg3In+L3xF7hr5kbmb4tjQt+mLHisR9UmeDAaRF5XncPLy9gurEJa8pVlhy0RU66FJbsT+GTdMY4kpdPArzqP9G7M6M4N8Tp/GI6tMvpqYzeDxWQUnWrc2+jDbzYAajcp+8kc/cKZo8R/9qgxMc3h3+DkFmO+4RpB4NYG3lkGJyse/4UMEzOWH+brrbH4e3swY2QE/VrWqaI3UgQHGfhgz6S7pirZcZ+ixaL54/AZZq09xvaY89TycueB7iFEdg+hrq8nZKcbNcGProKjK+B8jLFj7SbQdACE9oJG3cG3hP5YO/yQKxd7jd+SCye3XU7sKX8Z2+u2hZa3GY+g9sbfXwWTpNaaH3ae4q1fD3I+I4cHe4Ty3K0t8C1veWBJ0jYnSb4q2WuSuEp0zDlmrT3OqkNJuCrF4PB6PNQzlE4htS5XDEw5lpfwVxrJ35TXOvQLgeAeENzdWAa0uDxCw44/5MrEnuK/cAriNhv//kd+Ny6eu7gb37Ja3m5cPPcLtsqpDp2+yD9/2sf2mPN0CPbjX3eFE96gZvkP5CjfhJycJPmq5GB/5LEpl/hqcywLok9yMctM6yBfHuoZytD29a8cHpdrgtN/GrXw4zYby0vJxnOefnkJvzs89w7sOgXmq05kZx9yxbLVh7TFAskH8/598x4X4oznPP2MhN7yNuMblaf16sKkZ5t5b8URPt8Ug6+nG5NvC2Nkp0a4uFSwNLCDNHKcnST5quaAX1czcsz8tCuBeZtiOJyUhp+XO6O7NOKB7iE0rFXExTat4dzxy0n/5FY4e8R4zqLhrAWSLJCUC+fd4cXp8ODjRovYnlnjQ7osv39TFiTsvPyBeXIrZF0wnvOpe/lbUnB3o0vG1brlArTW/LI3kX8tPUDSxWzu69qIFweFUcu7WuUObE/fhG5gkuRFsbTWbDl+ji83x7D8QBJaawa2qsvYnqH0aOJfcgvv0lkjWS3/CnavBN9s8Cs0YMsrAOqFQ91wqBcBddsY/f32doNWZT6kr/6QcAHqe8GUCRBWx2itnzkEKUeNi9wAAS2vTOq1Qqv0w/DomXSmLtnPhqNnaVPfl38NC6djcIkzcZadtOTtgk2TvFJqMPA+4ArM0Vq/XdxrJcnbVkJqJl9vieXb7Sc5dymH4NpejOjYkBGdGhTdui9K5nk4vQ+S9uUt9xpJLrfQME3vOkZiK+pRI6jMd2TajMVi9JenJcLoQUAKBLpCHRfwdwHX/IStoFYIBLaCOmHQqJvx8Kp9XcI8dPoiH/1xjF/+TMDbw41/3NqS+7uH4FrRrpmiOFh3pbOyWZJXSrkCR4BbgHhgO3Cf1vpAUa+XJG8fsky5LNt3mu93nGTTsRS0hp5N/RnZuSGD2wRRvVo5b23PNcHZv+DMAWMET8EjFi7GG8MB87lWMy4u+tYHL3/jUb325Z+9al+5dPeqfCvYYgFzJpgyjYvN2emQfhrSThuJ/Ipl3iO/VQ5Gd8V5Dcm5kGy5/Ii5aJNvLXtOpjLzj6OsOJCEdzVX7u8RwqO9mxDg41E1J3TA7kpnY8sk3wOYqrUelLf+EoDW+q2iXi9J3v7En8/gx52nWLgjnrhzGfh4uDEkIoh7OjW8cmRORZlzjMJZqbFXfgCkJUFGivHIPA8U83fqWs24m9fV3Xi4uBv92a7VLv/s4m6s61wjiZsKJXRTJpizSo7Rs6bxDaNGvWuXDzwBf9rHhedtJ84x84+jrDuSjK+nG+N6NWZcr1D8vCrZ7y7snq0m8gZoAJwstB4PdKvicworaljLi6cHNOepm5ux7cQ5Fu6IZ8meBL7dfpLGAd7c06khwzo0oIFfBSeOcKsG/k2NR3EsuZCZanSR5Cf+wh8AuSZj0otcE1jMeUtT3vZCP7tUy2v9V897eBW/rFHPePjUK7k1/ly60V1hvqq74jrdsam1ZsPRs3y4+ijbTpzD37sakwaHcX/3YGqUdby7tMSdWlW35EcCg7TWj+StPwB01Vo/Veg144HxAMHBwZ1ii7qII+zKpWwzv+5NZOGOeLaeMOZmb9ewJoPC63FbeBCNA7xtHOF1ZoMkmWvRrD50hpl/HGXPyVTq+XryWN8m3NsluHzdadKn7hSku0ZUmbiUDH7Zm8iy/afZczIVgJZ1azA4vB6Dw+sRVq9G5bt0SnMDtUSPJKXxw854Fu9K4PTFLBrVrs7jfZsxolMDPNwqUAZYRsc4BVsmeTeMC68DgFMYF17HaK33F/V6SfKOLSE1k2X7TrNs/2m2x5xDawj19ypo4bdrWNP6Cf8GaIkmp2WzZE8Ci3bFs+/URVxdFP1aBDK8Y0MGtamLm2slRiPJOHenYOshlLcD72EMoZyrtS62s1KSvPNITstmxYEkftuXyOZjKZgtmqCantzULIBezQLo2dSfOr6elT+Rk7ZEs0y5rDiQxI8741n311lyLZq2DWoyvGMD7mxX33ojZZz03+9GIzdDCZu6kGFi5cEkVhxIYvPxFC5kGsMPm9XxoVdTf3o2C6B7Y39qepWzMBY4VUs0MyeX7THn+OXPRH7dm0hatpmgmp4M69CA4R0a0LxuFUz6cgN8E7oRSJIXdiPXojmQcJFNx86y8VgK20+cI9OUi4uC8AY16dnUaOV3CqmFt0cZBn85cEs0y5TLrrhUNh9PYfOxs+w+mYopV+NdzZXb2gYxvGMDujcu5a5ja7iBrmk4K0nywm7lmC3smruATUvWsalmKLsahGF2cUUpaBzgTXj9mrSp70ubvOU1tVYcqCWaY7bwZ3wqm4+lsOlYCjvjzpNttuCioG2DmvRoGkCPpv50Da1d/hvOxA1NkrywX1cl6Uvunmxr1pE9jzzL/oBQ9p+6QMKFyzcr1a/pSZsGlxN/6/q+1Fv6I66v2E9LVGtNclo2x89e4njyJY4np3M4KY3omPNkmnIBaB3kS4+m/vRs6k+XxrXLX8NdiEIkyQv7VYbulnOXcjiQcJH9CRfYl7c8cfZSQVe8m4uiXk1P6vtVp4Ffder7GT/nrwfV9Cz7jUFlZMq1cCnbTEJqFsfPphck8+NnL3Ei+RJp2ZdvgfVwc6FJoA9dQ2vRo6k/3Rr7V776oxCFSJIX9quCF04vZZs5dPoiBxPTSEjNzHtkcSo1k9MXs8i1XHnMGp5u1PauhqebK57uLni4ueLh7oKnuyue7q54uLng6e6Cp5srbq4uZOaYScs2k55lJj3bTNoVSxNZpmtja+BXnSaB3jQJ8KZxgDdNAn1oEuhN/ZrVq75fXdzQbFnWQIiSBQcX3ZIPLnkGJG8PNzqF1KZTyLUVHXMtmjNpWSSkZhV8AJxKzeRCpoksUy7ZZgtZplzSsswkp2WTk7eelbc05VrwquaGj4cbNTyNR4BPNUIDvAu2+XgYjzq+HjQJ8KFxgLf0owu7JEle2Na0aUVfOK1E7RdXF0VQzeoE1axOpxAr1U0XwkHZeeFu4fQiI42RMCEhRhdNSIhdjowRwlFJS17YXmSkJHUhqoi05IUQwolJkhdCCCcmSV4IIZyYJHkhhHBikuSFEMKJSZIXQggnJkleCCGcWKWSvFJqpFJqv1LKopTqfNVzLymljiqlDiulBlUuTCGEEBVR2Zuh9gHDgU8Kb1RKtQbuBdoA9YGVSqkWWuvcSp5PCCFEOVSqJa+1Pqi1PlzEU3cB32qts7XWJ4CjQNfKnEsIIUT5VVWffAPgZKH1+Lxt11BKjVdKRSulopOTk6soHCGEuDGV2l2jlFoJ1CviqSla68XF7VbEtiIL12utZwOzwagnX1o8Qgghyq7UJK+1HliB48YDjQqtNwQSKnAcIYQQlVBV3TVLgHuVUh5KqcZAc2BbFZ1LCCFEMSo7hPJupVQ80AP4RSn1O4DWej+wADgALAOekJE1Qghx/VVqCKXWehGwqJjnpgEVn95HCCFEpckdr0II4cQkyQshhBOTJC+EEE5MkryAqCgIDQUXF2MZFWXriIQQViITed/ooqJg/HjIyDDWY2ONdZDJtYVwAtKSv9FNmXI5wefLyDC2CyEcniT5G11cXPm2CyEciiT5G11wcPm2CyEciiT5G920aeDldeU2Ly9juxDC4UmSv9FFRsLs2RASAkoZy9mz5aKrEE5CRtcII6FLUhfCKUlLXgghnJgkeSGEcGKS5IUQwolJkhdCCCcmSV4IIZyY0tp+5s5WSiUDsZU4RABw1krh2JKzvA+Q92KPnOV9gLyXfCFa68CinrCrJF9ZSqlorXVnW8dRWc7yPkDeiz1ylvcB8l7KQrprhBDCiUmSF0IIJ+ZsSX62rQOwEmd5HyDvxR45y/sAeS+lcqo+eSGEEFdytpa8EEKIQiTJCyGEE3OqJK+U+pdS6k+l1G6l1HKlVH1bx1RRSqn/U0odyns/i5RSfraOqaKUUiOVUvuVUhallMMNd1NKDVZKHVZKHVVKTbZ1PBWllJqrlDqjlNpn61gqSynVSCn1h1LqYN7f1jO2jqkilFKeSqltSqk9ee/jdaufw5n65JVSvlrri3k/Pw201lpPsHFYFaKUuhVYrbU2K6X+A6C1nmTjsCpEKdUKsACfAP/QWkfbOKQyU0q5AkeAW4B4YDtwn9b6gE0DqwClVB8gHfhSax1u63gqQykVBARprXcqpWoAO4BhjvZ7UUopwFtrna6Ucgc2AM9orbdY6xxO1ZLPT/B5vAGH/QTTWi/XWpvzVrcADW0ZT2VorQ9qrQ/bOo4K6goc1Vof11rnAN8Cd9k4pgrRWq8Dztk6DmvQWidqrXfm/ZwGHAQa2Daq8tOG9LxV97yHVfOWUyV5AKXUNKXUSSASeNXW8VjJ34DfbB3EDaoBcLLQejwOmEycmVIqFOgAbLVxKBWilHJVSu0GzgArtNZWfR8Ol+SVUiuVUvuKeNwFoLWeorVuBEQBT9o22pKV9l7yXjMFMGO8H7tVlvfioFQR2xz2G6KzUUr5AD8Af7/qm7zD0Frnaq3bY3xb76qUsmpXmsNN/6e1HljGl34D/AK8VoXhVEpp70UpNRYYAgzQdn7xpBy/F0cTDzQqtN4QSLBRLKKQvD7sH4AorfWPto6nsrTWqUqpNcBgwGoXxx2uJV8SpVTzQqtDgUO2iqWylFKDgUnAUK11hq3juYFtB5orpRorpaoB9wJLbBzTDS/vguVnwEGt9X9tHU9FKaUC80fOKaWqAwOxct5yttE1PwAtMUZyxAITtNanbBtVxSiljgIeQErepi0OPFLobuBDIBBIBXZrrQfZNKhyUErdDrwHuAJztdbTbBtRxSil5gP9MEraJgGvaa0/s2lQFaSUuglYD+zF+P8O8LLW+lfbRVV+SqkIYB7G35YLsEBr/YZVz+FMSV4IIcSVnKq7RgghxJUkyQshhBOTJC+EEE5MkrwQQjgxSfJCCOHEJMkLIYQTkyQvhBBO7P8BYAELyEydIo4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"weights without regularization\")\n",
    "print(no_regularization_params)\n",
    "print(\"weights with regularization\")\n",
    "print(params)\n",
    "\n",
    "pyplot.scatter(x, y, c='r', label='true')\n",
    "pyplot.plot(xgrid, Xgrid_scaled@no_regularization_params, label='w/o regularization')\n",
    "pyplot.plot(xgrid, Xgrid_scaled@params, label='with regularization')\n",
    "pyplot.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using our helper script again to display both models with varying polynomial degree in an `ipywidget`, interact with the slider to explore the results. We set `regularized=True` this time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e0c720ccb75d4994b1bb51d650dd8ab5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=8, description='degree', max=15, min=1), Output()), _dom_classes=('widge…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interact_polyreg(max_degree, x, y, regularized=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to the regularization term, you won't see the wiggling curves even with a high degree polynomial, since the magnitude of the weights (the coefficients before each polynomial term) is much smaller."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ridge regression with scikit-learn\n",
    "\n",
    "You won't be surprised by now that **scikit-learn** offers a method to obtain a linear model with regularization: [`Ridge()`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html#sklearn.linear_model.Ridge). \n",
    "Look and ponder at the code cell below, where we call `Ridge()` and its `.fit()` method to fit the ridge regression model with our scaled matrix of features (also called the _design matrix_) and the $y$ data.\n",
    "The default value of the penalization parameter of $1.0$, so we don't specify it, but if you wanted to try a different value, you specify it as `alpha=value` in the argument list. Do explore the documentation page."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[12.48291164 28.26642615 11.09583654]\n",
      "-11.138315887456029\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAArgklEQVR4nO3deXxU1d348c83ISEEAiQhhDUJSMIWCUsC+LCjIK0oqMUKUXiqNm59bG1xK32pqPwen2pbatXSWLVYxwUX1LqyCCoKQlgEJEAAkxDAJATIQhaSmfP7YyYxSBJIZiaz5Pt+vfK6c++dued7JzPfe+bcc+8RYwxKKaX8U4CnA1BKKeU+muSVUsqPaZJXSik/pkleKaX8mCZ5pZTyY+08HUB93bp1M3FxcZ4OQymlfMrWrVuPG2OiGlrnVUk+Li6OjIwMT4ehlFI+RURyGlunzTVKKeXHNMkrpZQf0ySvlFJ+zKva5BtSXV1NXl4elZWVng5FtQEhISH06dOHoKAgT4eilEt4fZLPy8sjLCyMuLg4RMTT4Sg/ZoyhqKiIvLw8+vXr5+lwlHIJr2+uqaysJDIyUhO8cjsRITIyUn81qtZlsUBcHAQE2KcWi0s37/U1eUATvGo1+llTrcpigbQ0KC+3z+fk2OcBUlNdUoTX1+SVUspvLVr0Q4KvVV5uX+4imuRbKC4ujuPHj5+zvFOnTh6I5gcZGRncddddDa5rLGZ/sGzZMl566SVPh6FU8+TmNm95C/hEc01bZLVaCQwMbNZrampqSE5OJjk52U1RNVxmu3Yt/xgZYzDGEBDgXH3jtttuc+r1SnlETIy9iaah5S6iNfnzOH36NFdccQVJSUkkJiby+uuvn7W+oqKCGTNm8Nxzz53z2ieeeIKUlBSGDRvGQw89VLd89uzZjBo1iqFDh5Kenl63vFOnTjz44IOMGTOGjRs30qlTJxYtWkRSUhJjx44lPz//nDIefvhh0tLSmD59OvPnz2f9+vXMnDkTgKKiIqZPn86IESO49dZbqT8K2KOPPsqgQYOYNm0ac+fO5cknnwTg4MGDzJgxg1GjRjFhwgT27t173jILCwu59tprSUlJISUlhS+//BKAwsJCpk2bxsiRI7n11luJjY3l+PHjZGdnM3jwYO644w5GjhzJ4cOHG3yvGnvv77//foYMGcKwYcNYuHBhXUy1+7Bjxw7Gjh3LsGHDuPrqqzl58iQAkydP5r777mP06NEkJCTwxRdfNPm/V8rtliyB0NCzl4WG2pe7iE/V5Bf/51v2HC1x6TaH9OrMQ1cObXT9xx9/TK9evfjggw8AKC4urltXVlbG9ddfz/z585k/f/5Zr1u1ahVZWVls3rwZYwxXXXUVn3/+ORMnTuSFF14gIiKCiooKUlJSuPbaa4mMjOT06dMkJibyyCOPAPYkN3bsWJYsWcK9997Lc889xx/+8IdzYty6dSsbNmygQ4cOrF+/vm754sWLGT9+PA8++CAffPBB3QElIyODt956i+3bt1NTU8PIkSMZNWoUAGlpaSxbtoz4+Hi+/vpr7rjjDj799NMmy5w3bx53330348ePJzc3l8svv5zMzEwWL17M1KlTeeCBB/j444/POqDt27ePF198kWeffbbR96qwsPCc9/7EiROsXLmSvXv3IiKcOnXqnNjmz5/P3/72NyZNmsSDDz7I4sWLWbp0KWD/5bF582Y+/PBDFi9ezJo1axr93yvldrUnVxctsjfRxMTYE7yLTrqCjyV5T7j44otZuHAh9913HzNnzmTChAl162bNmsW9995LagP/kFWrVrFq1SpGjBgB2A8IWVlZTJw4kaeeeoqVK1cCcPjwYbKysoiMjCQwMJBrr722bhvBwcF1tfJRo0axevXqBmO86qqr6NChwznLP//8c95++20ArrjiCsLDwwHYsGEDs2bNqnvNlVdeWRfjV199xZw5c+q2UVVVdd4y16xZw549e+rWlZSUUFpayoYNG+r2c8aMGXXlA8TGxjJ27Ngm36sJEyac897X1NQQEhLCLbfcwhVXXFH3/tQqLi7m1KlTTJo0CYAFCxactT/XXHNN3fuZnZ3d4L4p1apSU12a1H/Mp5J8UzVud0lISGDr1q18+OGHPPDAA0yfPp0HH3wQgHHjxvHRRx8xb968c7reGWN44IEHuPXWW89avn79etasWcPGjRsJDQ1l8uTJdf2yQ0JCzmqHDwoKqttuYGAgNTU1DcbYsWPHRuNvqEtgY4O322w2unbtyo4dOxrdXkNl2mw2Nm7ceM6BpqlB4uu/vrH3Cmjwvd+8eTNr167ltdde4+mnn27wl0Zj2rdvDzT9firlT7RN/jyOHj1KaGgoN9xwAwsXLmTbtm116x555BEiIyO54447znnd5ZdfzgsvvEBZWRkAR44coaCggOLiYsLDwwkNDWXv3r1s2rTJbbFPnDgRi+PCio8++qiubXr8+PH85z//obKykrKysrrmkM6dO9OvXz/eeOMNwJ58v/nmm/OWM336dJ5++um6+dqDxPjx41mxYgVgr63Xlv9jjb1XDb33ZWVlFBcX89Of/pSlS5eec0Dq0qUL4eHhde3t//73v+tq9Uq1RT5Vk/eEXbt2cc899xAQEEBQUBB///vfz1q/dOlSbrrpJu69917++Mc/1i2fPn06mZmZXHLJJYD9pOrLL7/MjBkzWLZsGcOGDWPgwIF1TRbu8NBDDzF37lxGjhzJpEmTiHGcsU9JSeGqq64iKSmJ2NhYkpOT6dKlCwAWi4Xbb7+dxx57jOrqaq6//nqSkpKaLOepp57izjvvZNiwYdTU1DBx4kSWLVtWV/7rr7/OpEmT6NmzJ2FhYXXJvFZj79WBAwfOee9LS0uZNWsWlZWVGGP4y1/+ck48y5cv57bbbqO8vJz+/fvz4osvuuLtVMonSVM/qVtbcnKy+fGgIZmZmQwePNhDEfmvsrIyOnXqRHl5ORMnTiQ9PZ2RI0e6tIyqqioCAwNp164dGzdu5Pbbb7+gpiBP08+c8jUistUY02Dfaa3Jt1FpaWns2bOHyspKFixY4PIED5Cbm8t1112HzWYjODi4wW6mSin3ckmSF5GuwD+BRMAANwH7gNeBOCAbuM4Y03CjrGp1r7zyitvLiI+PZ/v27W4vRynVOFedeP0r8LExZhCQBGQC9wNrjTHxwFrHvFJKqVbkdJIXkc7AROB5AGPMGWPMKWAWsNzxtOXAbGfLUkop1TyuqMn3BwqBF0Vku4j8U0Q6AtHGmGMAjmn3hl4sImkikiEiGYWFhS4IRymlVC1XJPl2wEjg78aYEcBpmtE0Y4xJN8YkG2OSo6KiXBCOUkqpWq5I8nlAnjHma8f8m9iTfr6I9ARwTAtcUJZHnDp1imeffdbTYSilVLM5neSNMd8Dh0VkoGPRpcAe4D1ggWPZAuBdZ8vylMaSvNVq9UA0Sil14VzVu+Z/AIuI7ASGA/8PeByYJiJZwDTHvPu5YbzE+++/n4MHDzJ8+HBSUlKYMmUK8+bN4+KLLyY7O5vExMS65z755JM8/PDDwIXdtlcppdzJJf3kjTE7gIautrrUFdu/YG4aL/Hxxx9n9+7d7Nixg/Xr13PFFVewe/du+vXr1+SdDC/0tr1KKeUu/nXFa1PjJbrwVp6jR4+mX79+TT6nObftVUopd/GvJN8K4yXC2bfJbdeuHTabrW6+9rbBzbltr1JKuYt/3Wq4sXERnRwvMSwsjNLS0gbXRUdHU1BQQFFREVVVVbz//vtAy2/bq5RSruRfSd5N4yVGRkYybtw4EhMTueeee85aFxQUVDcu68yZMxk0aFDdOovFwvPPP09SUhJDhw7l3Xd9toORUspH+d+thi0Wt46XqPyf3mpY+Zq2dathN4+XqJRSvsS/mmuUUkqdxSeSvDc1KSn/pp815W+8PsmHhIRQVFSkXz7ldsYYioqKCAkJ8XQoSrmM17fJ9+nTh7y8PPQ2xKo1hISE0KdPH0+HoZTLeH2SDwoKOu/VpUoppRrm9c01SimlWk6TvFJK+TFN8kop5cc0ySullB/TJK+UUn7MJb1rRCQbKAWsQI0xJllEIoDXgTggG7jOGHPSFeUppZS6MK6syU8xxgyvd5Oc+4G1xph4YK1jXimlVCtyZ3PNLGC54/FyYLYby1JKKdUAVyV5A6wSka0i4hhUlWhjzDEAx7R7Qy8UkTQRyRCRDL2qVSmlXMtVV7yOM8YcFZHuwGoR2XuhLzTGpAPpYL+fvIviUUophYtq8saYo45pAbASGA3ki0hPAMe0wBVlKaWUunBOJ3kR6SgiYbWPgenAbuA9YIHjaQsAHftOKaVamSuaa6KBlSJSu71XjDEfi8gWYIWI3AzkAnNcUJZSSqlmcDrJG2MOAUkNLC8CLnV2+0oppVpOr3hVSik/pkleKaX8mCZ5pZTyY5rklVLKGRYLxMVBQIB9arF4OqKzeP3wf0op5bUsFkhLg/Jy+3xOjn0eIDXVc3HVozV5pZRqqUWLfkjwtcrL7cu9hCZ5pZRqqdzc5i33AE3ySinVUjExzVvuAZrklVKqpZYsgdDQs5eFhtqXewlN8kop1VKpqZCeDrGxIGKfpqd7zUlX0N41SinlnNRUr0rqP6Y1eaWU8mOa5JVSyo9pkldKKT+mSV4ppfyYy5K8iASKyHYRed8xHyEiq0UkyzENd1VZSinlT4rLq7HZ3DPEtStr8r8GMuvN3w+sNcbEA2sd80oppeqprLZy4wtf87s3vnHL9l2S5EWkD3AF8M96i2cByx2PlwOzXVGWUkr5C2MMi1buZmdeMT9J7OGWMlxVk18K3AvY6i2LNsYcA3BMu7uoLKWU8gv/+iqbt7bl8ZvL4pk+1EuTvIjMBAqMMVtb+Po0EckQkYzCwkJnw1FKKZ/w1cHjPPZBJtOGRHPX1Hi3leOKmvw44CoRyQZeA6aKyMtAvoj0BHBMCxp6sTEm3RiTbIxJjoqKckE4Sinl3fJOlvOrV7YTFxnKn69LIiBA3FaW00neGPOAMaaPMSYOuB741BhzA/AesMDxtAXAu86WpZRSvq7ijJVb/72V6hobz81PJiwkyK3lubOf/OPANBHJAqY55pVSqs0yxnD/2zvZc6yEv84dTv+oTm4v06U3KDPGrAfWOx4XAZe6cvtKKeXL/vnFd7y74ygLpycwdVB0q5SpV7wqpVQr2JB1nP/9KJOfJPbgzikDWq1cTfJKKeVmuUXl/OrVbcR3D+PJOUmIuO9E649pkldKKTcqP1ND2r8zsNkM6fNH0bF96w7joYOGKKWUmxhjuOfNnezPL+XFX4wmNrJjq8egNXmllHKTZZ8d4oOdx7h3xiAmJXjmOiBN8kop5QZrM/P54yd7mTmsJ7dO7O+xODTJK6WUi+3KK+ZXr2xnaK/O/PFnw1r1ROuPaZJXSikXOnyinJuWbyGiYzAvLEghNNizpz71xKtSSrlIcXk1v/jXFiqrrbxyyxi6dw7xdEhak1dKKVeoqrFy68sZ5BSd5h83jiI+OszTIQFak1dKKacZY7jvzZ1sOnSCpT8fzn9d1M3TIdXxj5q8xQJxcRAQYJ9aLJ6OSCnVhvxp1X7ecdyTZvaI3p4O5yy+X5O3WCAtDcrL7fM5OfZ5gNRUz8WllGoTXtucy9PrDnB9St9WvSfNhfL9mvyiRT8k+Frl5fblSinlRp/tL2TRO7uZmBDFo7MTPdpVsjG+n+Rzc5u3XCmlXODbo8Xc8fJWBkaH8WzqSIICvTOdemdUzRET07zlSinlpKOnKrjpX1vo3CGIF3+RQqdWvulYc7hiIO8QEdksIt+IyLcistixPEJEVotIlmMa7ny4DViyBEJDz14WGmpfrpRSLlZSWc1N/9pCeZWVF3+RQrQX9IVviitq8lXAVGNMEjAcmCEiY4H7gbXGmHhgrWPe9VJTIT0dYmNBxD5NT9eTrkopl6ustnL7y1s5UFDGshtHMahHZ0+HdF5O/8YwxhigzDEb5PgzwCxgsmP5cuzDAt7nbHkNSk3VpK6Ucqtqq41fvbKNLw8U8efrkhg3wHv6wjfFJW3yIhIoIjuAAmC1MeZrINoYcwzAMe3eyGvTRCRDRDIKCwtdEY5SSrmU1Wa4+/UdrMks4NHZiVwzso+nQ7pgLknyxhirMWY40AcYLSKJzXhtujEm2RiTHBXlmfstK6VUY2w2wwNv7+T9ncd44CeDuHFsrKdDahaX9q4xxpzC3iwzA8gXkZ4AjmmBK8tSSil3M8bwyPt7WJGRx12XxnPrpIs8HVKzuaJ3TZSIdHU87gBcBuwF3gMWOJ62AHjX2bKUUqo1PblqH//6Kptbxvfj7sviPR1Oi7iic2dPYLmIBGI/aKwwxrwvIhuBFSJyM5ALzHFBWUop1SqeWXeAZ9YdZO7oGBZdMdgrr2a9EK7oXbMTGNHA8iLgUme3r5RSre3FL7/jiU/2MXt4Lx7z0tsVXCjfv+JVKaVcaMWWwyz+zx4uHxrNk3OSCAzw3QQPmuSVUqrOe98c5b63dzIxIYqn5o6gnZfej6Y5fH8PlFLKBVbvyee3r+8gJTaCf9wwivbtAj0dkktokldKtXnr9hVwp2UbQ3t15vn/TqZDsH8keNAkr5Rq4z7efYy0lzKIj+7E8ptGExYS5OmQXEqTvFKqzXpn+xHufGU7F/fuwiu/HEvX0GBPh+Ry3nsTZKWUcqNXN+fy+5W7GNMvgucXpNDRi+8J7wz/3CullGrCCxu+45H39zB5YBTLbhhFSJD/tMH/mCZ5pVSb8sy6AzzxyT5mDO3BX+cO95teNI3RJK+UahOMMTy5ah/PrDvI7OG9eHJOkl/0gz8fTfJKKb9XezfJF7/MZu7ovjw2+2Kfv5L1QmmSV0r5NavN8Id3dvHq5sP8YlwcD84c4tP3omku//+topTybxYLxMVBQIB9arHUraqx2vjdih28uvkwd065qM0leNCavFLKl1kskJYG5eX2+Zwc+zxQed31/Pq17XzybT73XD6QO6cM8GCgnqM1eaWU71q06IcEX6u8nOOP/C/Xp29i1Z58Hr5ySNMJvolfAv5Aa/JKKd+Vm3vOogMRffjFlLsp/L6Ev6eOYkZij8Zf38QvAVJT3RBw63PF8H99RWSdiGSKyLci8mvH8ggRWS0iWY5puPPhKqVUPTExZ81+3Wco197wBBUhobyWdknTCR4a/SXAokUuDtRzXNFcUwP8zhgzGBgL3CkiQ4D7gbXGmHhgrWNeKaVcZ8kSCA0F4J0hk7nx54/RraKYlUkwvG/X87++gV8CTS73Qa4Y/u8YcMzxuFREMoHewCxgsuNpy4H1wH3OlqeUUnVSUzEGnn7lC/407ErG5u/nH1Oi6TJ/7oW9PibG3kTT0HI/4dI2eRGJwz7e69dAtOMAgDHmmIh0b+Q1aUAaQIwfvbFKKferttr4fftE3hjWlWtG9Obxa39CcLtmNFAsWXJ2mzzYfxksWeL6YD3EZb1rRKQT8BbwG2NMyYW+zhiTboxJNsYkR0VFuSocpZSfK66o5r9f3MwbW/P49aXx/Om6pOYleLCfXE1Ph9hYELFP09P95qQruKgmLyJB2BO8xRjztmNxvoj0dNTiewIFrihLKaXyTpZz07+2cKjwNE/OSeJno/q0fGOpqX6V1H/MFb1rBHgeyDTG/LneqveABY7HC4B3nS1LKaW25Z7k6me/4lhxJS/dNNq5BN8GuKImPw64EdglIjscy34PPA6sEJGbgVxgjgvKUkq1UcYY/r0ph0ff30OPLiFYbhlDQnSYp8Pyeq7oXbMBaOxmEJc6u32llCo/U8OilbtZuf0IUwZGsfTnI+gS6l9jsbqLX9zW4HRVDXe9up1vjxZ7OhSllItlHz/NNc9+xTs7jvDbaQk8vyBFE3wz+MVtDfZ+X8pn+wt575ujzBrei99NG0hMZKinw1JKOWn1nnx+u2IHgQHCv34xmkkJ2gOvucQY4+kY6iQnJ5uMjIwWvba4opp/fHaQF778DqvNkDomll9NHUC3Tu1dHKVSyt2sNsOfVu3j2fUHubh3F55NHUnfCK24NUZEthpjkhtc5y9JvlZ+SSVL12SxIuMwIe0C+OXE/twyoT+d/HQkdqX8TVFZFXe9tp0vDxQxd3RfHrpyqF8PtO0KbSrJ1zpYWMafVu3jw13fE9kxmP+ZOoB5Y2Kbf7GEUqrVbM89yR2WbRSdPsNjsxK5LqWvp0PyCW0yydfacfgU//fRXjYeKqJvRAcWTh/IlcN6EdBGxndUyhfU7x4Z3TmEZTeMIrF3F0+H5TPadJIH+wfo86zjPP7RXjKPlTCoRxh3T0tg+pDoNjcUmFLepqCkknvf2sn6fYVMHhjF0p8Pp2tosKfD8iltPsnXstkM/9l5lKVrsvju+GkSe3fm7ssSmDqouyZ7pTzgg53HWPTOLiqrrfz+p4O5cWysfhdbQJP8j9RYbbyz4yhPrc0i90Q5SX26cPe0BCYlROkHTKlWUFxRzcPvfcvK7UdI6tOFP/98OBdFdfJ0WD5Lk3wjqq023t6Wx1NrD3DkVAUjY7ry22kDGTcgUpO9Um7y1YHjLHzjG/JLq/ifqQO4c8oAggK1Q4QzNMmfx5kaG29sPczTnx7gWHElo+MiuHtaApdcFNnqsSjlryqrrfzx43288OV39O/Wkb/8fDhJFzJ6kzovTfIXqKrGyutbDvPMugPkl1Qxtn8Ed02N55KLtGavlDN2HynmN6/v4EBBGQsuieX+nwymQ7D2fXcVTfLNVFlt5dXNuTy7/iCFpVUM79uVO6cM4NJB3bXrpVLNUGO1seyzgyxdk0Vkp2Ce+FkSE/XWBC6nSb6FKqutvLk1j2WfHSTvZAUDo8O4ffJFzBzWk3bahqhUkzKyT/CHd3az9/tSrkzqxaOzhmrXSDdpKslrpmpCSFAgN4yNZf3Cyfzl50nYjOE3r+9g6p8+w/J1DpXVVk+HqPyBxQJxcRAQYJ9aLJ6OyCknTp/h3je/4WfLNlJSUc2yG0bxt7kjNMF7iNbkm8FmM6zOzOfZdQf4Jq+Y7mHt+eWE/swbE0NHZ+6NY7HAokWQm2sfJX7JEr8ejkzVY7E0PJC0D44zarMZVmQc5vGP91JWWcPNE/px19R4574b6oK4vblGRF4AZgIFxphEx7II4HUgDsgGrjPGnGxqO96e5GsZY/jqYBHPrDvAVweL6BoaxPxL4rhxbCxRYc2866UffclVC8TFQU7OuctjYyE7u7WjabFvjxbzh3d2sz33FKP7RfDY7EQdtakVtUaSnwiUAS/VS/J/BE4YYx4XkfuBcGPMfU1tx1eSfH3bck/y7LqDrMnMJzgwgFnDe3HzhH4M6tH5wjbgJ19y1UIBAdDQd1AEbLbWj6eZSiur+fPq/Sz/Kpvw0GB+/9PBXDOyt/ZGa2WtcuJVROKA9+sl+X3AZGPMMRHpCaw3xgxsahu+mORrHSos48Uvs3lzax4V1VbGDYjk5vH9mJxwnh45Pv4lV07y0YO8MYb3dx7j0ff3UFhWReqYGO6ZPkhHbPIQTyX5U8aYrvXWnzTGhDfwujQgDSAmJmZUTkMfeB9yqvwMr24+zPKvsvm+pJL+UR25aVw/rh3Zp+F+wT76JVcu4oPNdfXv7Hpx7y48NjtRL2ryMK9O8vX5ck3+x6qtNj7cdYznN3zHzrxiuoYGMW90DAv+K47oziE/PNEHv+TKxXzkxHtWfilPrtrHJ9/mE9kxmN9cFs+8MbEE6rUjHqfNNR5kjCEj5yTPf/Edq/Z8T4AIlyf2IHVMDJf0d1xJ6yNfctU25Z0sZ+maLN7elkdocDt+OaE/N0/op6OteZGmkrw7/0vvAQuAxx3Td91YltcSEVLiIkiJiyC3qJyXNmbz5rY8Pth5jP7dOjJvTAzXzp5DuCZ15WWOl1XxzLoDWDblgsBN4/pxx5QBRHTU/u6+xFW9a14FJgPdgHzgIeAdYAUQA+QCc4wxJ5rajsdq8q1ck66stvLR7mNYNuWSkXOS4HYBzLy4J/PGxDAqNlx7JiiPKq2s5rkvvuP5Lw5RUW1lzqi+/PqyeHp17eDp0FQj9LYGTfFwm/je70t45etcVm47QmlVDQOjw0gdG8PsEb3pHHKBPRW0uUe5QMUZK5avc3hm3QFOllfz04t78NtpAxnQXe/z7u00yTfFS3q3lJ+p4T/fHMXydS4784rpEBTIzGE9uWZkH8b0i2i8G6aeuFVOKiyt4qWN2by8KYeT5dVMiO/GPZcPZFifrp4OTV0gTfJN8cJ+6jvzTvHK17m8v/MYZVU19O7agatH9Obqkb3PHT3HSw5SyvccKCjln198x9vbj1BttXHZ4GjSJvYnJS7C06GpZtIk3xQvTpIVZ6ys2vM9b287whdZhdgMJPXtyrUjezNzWC/7CTAvPEgp72WMYdOhEzz3xSE+3VtA+3YB/GxUH24e34/+Ovyez9Ik3xQfae4oKKnkvW+O8ta2I2QeKyEoUJgysDvXpD/KlE0f0t5ac/YLvOAgpbxH7XUbz31xiN1HSojsGMz8S+K4YWwMkZ2aeb8lV9NzSk7TJH8+PvYh23O0hJXb83hnx1EKS6voUlnG9P0b+cm+LxmXs4P27YO97iCl3KiJz29BaSVvbzvCvzfmcORUBf2jOnLL+P5cM7I3IUFeMDKTj1SyvJ0meT9VY7Xx5cEi3nl7A2sKrJQGdyDsTAVTowL5yZVjmZTQXYdY83cNJMmqTp359P+e443Qfny2vxCrzTC6XwRpE/oz1dtGN/Pi5lJfokm+DaiqsfLVgSI+2n2M1XvyOVleTYegQKYMimJGYk+mDuquVyj6I0eSNMC30RfxZuKlvDN0Mqc6dKZH5xCuGdmba0f1OfeEvbfQc0ou4akrXtsOL2juad8ukCmDujNlUHdqrDa+/u4EH+0+xiff5vPhru8JbhfAxPhuTB/Sg4kJUfToEnL+jSqvd/x4Me8kz+LNiy9jb/d+BNecYXrWJubsWsP4Axnef1+ZmJiGa/IxMa0fi5/SmryzvLxN0WozbMs9yUe7vufj3cc4WlwJwKAeYUxMiGJSQhTJceG0b+dEs44XHOSc4mPxnzx9hnX7Cvho9/es232UmoBAko7uY86uNVyZ+Tldqk77TnOHl39/fIU217iTD7UpGmPY+30pn+8v5LP9hWzJPkG11dAhKJD/uiiSSQOjmBgfRVy3jhe+UV//kvpI/IcKy1ibWcDqzHwysk9gMxDduT2zg4v52dIHiD+S9cOTvTD+JvnYQdYbaZJ3Jx9uUzxdVcOmQ0V85kj6OUX2RBcbGcrE+CjG9I8gOTai6aYdHzrINchL46/9BbZmTz6rM/M5VHgagME9OzNtcHcuGxJNYq8u9pOonk6Sni5faZJ3Ky9NEi2Rffw0n2cV8tm+QjYeKqL8jBWAvhEdSI6NIDkunJS4CAZEdfqhh4YPH+QAr4r/WHEFW7JP8tm+Qj7daz95HhQojO0fybQh0Uwd1J0+4aGtGtN5+cgvIX+nSd6d/PRDXm21kXmshC3ZJ8nIPsGW7JMcL6sCoEuHIJJjw0mOiyB5YRoXf7OBkJozZ2/AVw5yHjpI22yGrIIytmSfqHt/j5yqAOzv79RB3blscDQTE7oRdqE3qvMEP6rk+DJN8u7WBn6uGmPIKSonI6c26Z/goKMJIdBmpf+JPAYXfGf/Kz7KkPt+RdSCud5/22RXHKQv4P9fWW1l15FiR1K3v4cllfarlKPC2pMSF05yrH3cgcE9w2gXGOCqPXQvL/ol1JZpklduUVRWxdack+z6+Esyd+wns1MPjnTpXrc+smMwg3t2ZnDPMIb06sygHp2Ji+zofRdoOXOQ/tFBojogkJye/cn6/WPsTxhBVkEpWfllHDpeRrXV/l0b0L3TWUm9b0QH7z8YNkZr8l7Bo0leRGYAfwUCgX8aYx5v7Lma5H1fcXk1e46VkFn7930J+/PLOFPzQ62uW6f2xER0ICYilJiIUPo6pjGRoUSHhXjXFZkNsNkMJ8vPkF9SRfZ1C8gyIezvFsOBbjEciuhNdaC9eUUE+oaHkhDdifjoMEbFhDMqNpxwfxpZyU+bK32Nx5K8iAQC+4FpQB6wBZhrjNnT0PM1yfunaquNQ4Wn2ZdfyuET5eQWlZN7wv53rLgCW72PYHBgAH3CO9CjSwjhHYOJCA0mPDTI/rhjMOGh9mnX0CAiOgbTISjQ6VqwzWaorLFSccZKRbWV01VWCkoryS+pIr+kkoISx+PSSgpKqigorayrlQOIsdH3VD4Jx3MZUJRLwvFcEooOc1HOXu/71eIObaC50tt58orX0cABY8whRyCvAbOABpO88k9BgQEM7BHGwB5h56w7U2Pj6KkKDp/8IfEfPlFOQUkVmcdKOHn6DKcqqhts9gX7QSG4XQBBgUK7wACCAhzTQCEoMIB2gUK7gACCAwOwGkPFGSuV1fZkXlFtT+xVNU23HXcOaUd05xCiO4cwpn9H++Ow9kR3DqHP/OsYsHsLHWqqzn5RbCy0hQQP9oSuSd1ruTvJ9wYO15vPA8a4uUzlQ4LbBRDXrWOTF2BZbYbiimpOlp/h5OkznDh9hpPlZzhxuppTFWeorjFUW23U2GxUWw01VhvVNkN1jY0am31dtdVGcEAg4aHBdAgOpENQAB2CAgkJDqRDkOMvOJCQoEBCgwPpHhZCdOf2dA8Labo2vvA2SNsG9e/0HBpqr836Cq2J+zV3J/mGfkefVScTkTQgDSBG71ehGhAYIEQ4mmuI8nQ0P1KbDH01Sf64TT0nxz4PvrMPqknu7qeVB/StN98HOFr/CcaYdGNMsjEmOSrK277ByidYLPZeHgEB9qnF0rrlp6bae5LYbPapLyXHRYvOPmkK9vlFizwTj3I5d9fktwDxItIPOAJcD8xzc5mqLdGaqHNyc5u3XPkct9bkjTE1wK+AT4BMYIUx5lt3lqnaGK2JOqexJlJtOvUbbr+szhjzoTEmwRhzkTHGh85GKZ+gNVHnLFliP1Fcn6+dOFZN8pFrp5VqhNZEnZOaar9wKTbWfvVWbKxeyORnNMkrz3PmxKnWRJ3nyyeO1XlpkleeVXviNCfHfqOr2hOnF5rotSaqVJP0BmXKs/QGV0o5ranbGmhNXnmWnjhVyq00ySvP0hOnSrmVJnnlWXriVCm30iSvPEtPnCrlVu6+rYFS56e3qlXKbbQmr5RSfkyTvFJK+TFN8kop5cc0ySullB/TJK+UUn5Mk7xSSvkxTfJKKeXHnEryIjJHRL4VEZuIJP9o3QMickBE9onI5c6FqZRSqiWcvRhqN3AN8I/6C0VkCPbxXIcCvYA1IpJgjLE6WZ5SSqlmcKomb4zJNMbsa2DVLOA1Y0yVMeY74AAw2pmylFJKNZ+72uR7A4frzec5lp1DRNJEJENEMgoLC90UjlJKtU3nba4RkTVAjwZWLTLGvNvYyxpY1uDoJMaYdCAd7IOGnC8epZRSF+68Sd4Yc1kLtpsH9K033wc42oLtKKWUcoK7mmveA64XkfYi0g+IBza7qSyllFKNcLYL5dUikgdcAnwgIp8AGGO+BVYAe4CPgTu1Z41SSrU+p7pQGmNWAisbWbcE0OF9lFLKg/SKV6WU8mOa5JVSyo9pkldKKT+mSV6BxQJxcRAQYJ9aLJ6OSCnlIjqQd1tnsUBaGpSX2+dzcuzzoINrK+UHtCbf1i1a9EOCr1Vebl+ulPJ5muTbutzc5i1XSvkUTfJtXUxM85YrpXyKJvm2bskSCA09e1loqH25UsrnaZJv61JTIT0dYmNBxD5NT9eTrkr5Ce1do+wJXZO6Un5Ja/JKKeXHNMkrpZQf0ySvlFJ+TJO8Ukr5MU3ySinlx8QY7xk7W0QKgRwnNtENOO6icDzJX/YDdF+8kb/sB+i+1Io1xkQ1tMKrkryzRCTDGJPs6Tic5S/7Abov3shf9gN0Xy6ENtcopZQf0ySvlFJ+zN+SfLqnA3ARf9kP0H3xRv6yH6D7cl5+1SavlFLqbP5Wk1dKKVWPJnmllPJjfpXkReRREdkpIjtEZJWI9PJ0TC0lIk+IyF7H/qwUka6ejqmlRGSOiHwrIjYR8bnubiIyQ0T2icgBEbnf0/G0lIi8ICIFIrLb07E4S0T6isg6Ecl0fLZ+7emYWkJEQkRks4h849iPxS4vw5/a5EWkszGmxPH4LmCIMeY2D4fVIiIyHfjUGFMjIv8HYIy5z8NhtYiIDAZswD+AhcaYDA+HdMFEJBDYD0wD8oAtwFxjzB6PBtYCIjIRKANeMsYkejoeZ4hIT6CnMWabiIQBW4HZvvZ/EREBOhpjykQkCNgA/NoYs8lVZfhVTb42wTt0BHz2CGaMWWWMqXHMbgL6eDIeZxhjMo0x+zwdRwuNBg4YYw4ZY84ArwGzPBxTixhjPgdOeDoOVzDGHDPGbHM8LgUygd6ejar5jF2ZYzbI8efSvOVXSR5ARJaIyGEgFXjQ0/G4yE3AR54Ooo3qDRyuN5+HDyYTfyYiccAI4GsPh9IiIhIoIjuAAmC1Mcal++FzSV5E1ojI7gb+ZgEYYxYZY/oCFuBXno22aefbF8dzFgE12PfHa13IvvgoaWCZz/5C9Dci0gl4C/jNj37J+wxjjNUYMxz7r/XRIuLSpjSfG/7PGHPZBT71FeAD4CE3huOU8+2LiCwAZgKXGi8/edKM/4uvyQP61pvvAxz1UCyqHkcb9luAxRjztqfjcZYx5pSIrAdmAC47Oe5zNfmmiEh8vdmrgL2eisVZIjIDuA+4yhhT7ul42rAtQLyI9BORYOB64D0Px9TmOU5YPg9kGmP+7Ol4WkpEomp7zolIB+AyXJy3/K13zVvAQOw9OXKA24wxRzwbVcuIyAGgPVDkWLTJh3sKXQ38DYgCTgE7jDGXezSoZhCRnwJLgUDgBWPMEs9G1DIi8iowGfstbfOBh4wxz3s0qBYSkfHAF8Au7N93gN8bYz70XFTNJyLDgOXYP1sBwApjzCMuLcOfkrxSSqmz+VVzjVJKqbNpkldKKT+mSV4ppfyYJnmllPJjmuSVUsqPaZJXSik/pkleKaX82P8Hv6QEveVo2RgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "\n",
    "model = Ridge().fit(X_scaled[:, 1:], y) # Ridge() by default fits an intercept\n",
    "y_pred_sklearn = model.predict(Xgrid_scaled[:,1:])\n",
    "\n",
    "print(model.coef_)\n",
    "print(model.intercept_)\n",
    "\n",
    "pyplot.scatter(x, y, c='r', label='true')\n",
    "pyplot.plot(xgrid, y_pred_sklearn, label='sklearn ridge regression')\n",
    "pyplot.legend();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `model` variable we created above using `Ridge()` stores the model weights and the intercept in the `.coef_` and `intercept_` attributes. \n",
    "Compare the values with those we obtained above using the regularization and 3rd-order polynomial. Pretty close!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we've learned\n",
    "\n",
    "- Polynomial regression is a special case of multiple linear regression.\n",
    "- Pick too low-order a polynomial (e.g., a line) and the model could _underfit_ the data, giving large residuals (training error).\n",
    "- Pick too high a polynomial order, and you will get _overfitting_, as the model tries to fit the noise in the data resulting in wiggles.\n",
    "- Regularization is used to reduce or avoid overfitting.\n",
    "- Adding the sum-of-squares of the weights to the loss function _penalizes_ large weights and controls overfitting. \n",
    "- Linear regression with this so-called Tikhonov regularization is called _ridge regression_.\n",
    "- With **scikit-learn**, you can use `Ridge()` to fit a linear model with regularization. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='https://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}