notebooks_en/1_Linear_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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear regression by gradient descent\n",
    "\n",
    "This module of _Engineering Computations_ takes a step-by-step approach to introduce you to the essential ideas of deep learning, an algorithmic technology that is taking the world by storm. \n",
    "It is at the core of the artificial intelligence boom, and we think every scientist and engineer should understand the basics, at least. \n",
    "\n",
    "Another term for deep learning is deep neural networks. \n",
    "In this module, you will learn how neural-network models are built, computationally. \n",
    "The inspiration for deep learning may have been how the brain works, but in practice what we have is a method to build models, using mostly linear algebra and a little bit of calculus. \n",
    "These models are not magical, or even \"intelligent\"—they are just about _optimization_, which every engineer knows about!\n",
    "\n",
    "In this lesson, we take the first step of model-building: linear regression. \n",
    "The very first module of the _Engineering Computations_ series discusses [linear regression with real data](http://go.gwu.edu/engcomp1lesson5), and there we found the model parameters (slope and $y$-intercept) analytically. \n",
    "Let's forget about that for this lesson. \n",
    "The key concept we introduce here will be _gradient descent_. Start your ride here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\n",
    "\n",
    "This lesson is partly based on a tutorial at the 2019 SciPy Conference by Eric Ma [1]. He begins his tutorial by presenting the idea of _gradient descent_ with a simple quadratic function: the question is how do we find this function's minimum?\n",
    "\n",
    "$$f(w) = w^2 +3w -5$$\n",
    "\n",
    "We know from calculus that at the minimum, the derivative of the function is zero (the tangent to the function curve is horizontal), and the second derivative is positive (the curve slants _up_ on each side of the minimum). \n",
    "The analytical derivative of the function above is $f^\\prime(w) = 2w + 3$ and the second derivative is $f^{\\prime\\prime}(w)=2>0$. Thus, we make $2w+3=0$ to find the minimum.\n",
    "\n",
    "Let's play with this function using SymPy. We'll later use NumPy, and make plots with Matplotlib, so we load all the libraries in one place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "import numpy\n",
    "\n",
    "from matplotlib import pyplot\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run this SymPy method to get beautiful typeset symbols and equations (in the Jupyter notebook, it will use [MathJax](https://en.wikipedia.org/wiki/MathJax) by default): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll define the Python variable `w` to be a SymPy symbol, and create the expression `f` to match the mathematical function above, and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAG8AAAAVCAYAAABIfLDHAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADoUlEQVRoBe2Z61HbQBCAhYcCCOnA6QBSAnQQkgoCHYThF/7HOB0AFSShA0IFAXdgOgjjDsj3nXUaSSNsZMuyTbwz691b3WNft3eSt56fn5N5oNfr7TD+LJ2jm9KvyEfzzLsZO90D29O7TO3RJ1AnsRf8JfwD+CHKNnQxHug0MO0xATvIzdOH7yLby8k27AI80ETw3HX3C9BtM+U0D3jmNYnn5+d9cNjknP/rXPhxr2w7sh2wq7yJMy/Lj7RUfkKwnwlbZtDBS1M8g71M2fZc/t2yKk0sd4fe2jBIJ5MXgn8bC17qNM+7ffil3DRTQ0+hMXgJvMl0Cz0Cb7R8jeAp1dX7wyOo/hfYEfzbSPCYzOzWaYfQJG1LXbAWMMbLjxeeq1oDx52PIV6gDFYMVNxxvs5E2QxTL2XIADuOXlp57gsLkxs4Xw8u4fdE+FMwZg1sLbA0xPJQayCdLS9mZchMB6NPxtt+S9DEzvOdTmdLM8BpWenKhAtmWNNd9i6/DDLLpmCCvSkIwcNAnW/JeQ8OaWclC96ddAINwYDqjOxso11wFs9WBtDNEqyu6p+3aWZ72zYOvY2L+hobq5xnXrjAxLJ5huA7D/6AGpsHzwoHRRjB2P4YBatGscXy/Q29TDgNvS/puC72GrSfxgb0KBIf4E3KpAPjzjJogheO8lllx1sfCvS3NOmQ2pcRx7cB6OhBr8Ee9j9ADQ7lE7o29qLrIehmCQCvz/V/OAIsm48IwzaE/wxegAGQa6jRjze28QMyOZ0otmtTxqtAyKDS4F3bPK86MyfevkrzhCbz3IA64BfUEr8Qe5lbP92B0teCry/R968dYwAPGNfd5kfDEqiZ6cLZ2QCvc0c8Ky8wRD4XMGdVcKIevipYxmsBY0w25yjra9nUFo0OrwvQRu1lPv3Y2McJ5rPa7UJfmnOnQ4cIlkwzOwQzFSor7Dqe66Cyc9LuSyfeeC2RJuE0WHV7vVNU2REr0yAfvC6dy+dYlewLzikEdJqXWnxu4sUymV9WRwh5vatsq5Ity94r/Fz1t5oVJNjRCSaNf9xNWaQZaCejrEEBUll2eYnyFaLexgr6oXMsj74uGNwIq26vHz0K76a0vUEL4avLll+nI+Q6a+Rf0PPvGnwClRXeAWk3Duigs2c681SG8SZd/pOSyVf5YXoV7FXnlwD91N2EFNxIxsHPkCEJC8Gzx7IBxeYK3rL1b3P9fNlsc91Ja5lVIbMmddo8S5J/Uya9Zk7UnmQAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle w^{2} + 3 w - 5$"
      ],
      "text/plain": [
       " 2          \n",
       "w  + 3⋅w - 5"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w = sympy.Symbol('w', real=True)\n",
    "\n",
    "f = w**2 + 3*w - 5\n",
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAawAAAEeCAYAAADWyiHSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1+ElEQVR4nO3deVxU9eL/8dewqqi4ASKLgiAiihuupaWImZVWltnVotWWe1tuZWX39stWzbK9a1natWtl2aJlLqVmqbnhrrigiAIiO6LsMOf3hze+19xQYc4M834+Hj3SMzOctweY95xzPudzLIZhICIiYu9czA4gIiJSEyosERFxCCosERFxCCosERFxCCosERFxCCosERFxCCosERFxCCosERFxCCosERFxCG5mBxC5UBaLpQFwLTAAaAOUADuBHw3D2GVmNhGpOxZNzSSOxGKxTAKuA1YCm4AsoAHQARj03z8/bhjGdpMiikgdUWGJQ7FYLNcYhvHjOR73BYINw0iwYSwRsQEVljgki8USahhGstk5RMR2VFjikCwWy29AALAR+A1YZRjGDnNTiUhdUmGJw7JYLB5AL+BK4D6gsWEYLUwNJSJ1RqMExSFZLJbLOTlKcADQDFgIrDIzk4jUrfPtYWn3S+ySq6srMTExTJw4keHDh+Ph4WF2JACGDRvGkiVLzI4h4kgsNX2i9rDEIeXm5rJmzRp+++033nnnHVxcXOjXrx8vvviiqblycnJMXb9IfabCEofUrFkzQkNDSU1NJS0tjd9//52KigqzY4lIHVJhiUNq3749ERERDBgwgPvvv59PPvnEbg4LikjdUGGJQ0pKSsLFRVNhijgT/caLQ3nppZfIy8s7a1mtWLGChQsX2jiViNiC9rDEoXTp0oXrrruOBg0a0KNHD3x8fCgtLSUpKYmtW7cyZMgQnnnmGbNjikgdUGGJQ/n6669Zs2YNU6dOxdfXl4yMDJo2bcq4ceOYMWMGDRs2NDuiiNQRFZY4lE2bNnHo0CE+++wzfvnll1MeKykpUWGJ1GMqLHEo999/P8OGDSM5OZmYmJjq5YZhYLFYSE7WfLgi9ZVmuhCH9MADDzB9+nSzY5wmJiaGhATd2UTqr9KKKhq4u9bml6zxTBcaJSgOyR7LSsQZ/P3Lrdz9742mrFuFJSIiNZKcfYIlu47S0b+JKetXYYmISI3M+C0Zd1cX7ugfYsr6VVgiInJemYWlfLs5ndExgfg08TQlgwpLRETOa9bqg1RarYwf0N60DCoscRp33XUXvr6+dO7cuXrZhAkT6NixI9HR0dxwww0UFBRUPzZ58mTCwsKIiIhg6dKlJiQWsQ/HSir4bP1hroluQ3DLRqblUGGJ07jjjjtOu7liXFwcO3fuZPv27XTo0IHJkycDkJiYyNy5c9m1axdLlizhwQcfpKqqyozYIqabs+4QJ8oquW9gqKk5VFjiNAYOHEiLFi1OWTZ06FDc3E5eP9+3b1/S0tIAWLBgAWPGjMHT05OQkBDCwsLYsGGDzTOLmK20oopP1qQwsIMPnQO8Tc2iwhL5r1mzZnH11VcDkJ6eTlBQUPVjgYGBpKenmxVNxDRfb0oj50QZ919h7t4VaGomEQBefvll3NzcGDt2LHByqqc/s1jOfEH+jBkzmDFjBgDZ2dl1F1LExqqsBh+tSqZrUDP6hbY0O472sERmz57NwoUL+eyzz6pLKTAwkNTU1OrnpKWl0aZNmzO+fvz48SQkJJCQkICPj49NMovYwuKdGRzKLeaBK0LP+oHNllRY4tSWLFnCq6++yvfff0+jRv83+mnEiBHMnTuXsrIyDh48SFJSEr179zYxqYhtGYbBD9uOEOrjxdBOrc2OA+iQoDiRW2+9lZUrV5KTk0NgYCDPP/88kydPpqysjLi4OODkwIsPPviAqKgoRo8eTadOnXBzc+P999/H1bVWJ/wUsWsr92azdFcmb4yOxsXF/L0r0GztIrVKs7VLfWAYBjd9sJaMghJWThiEh1udHozTbO0iInJx1h/MY9OhfO67on1dl9UFsZ8kIiJiF97/ZT+tGntyS6+g8z/ZhlRYIiJSbVtqAauScrhnQEht36jxkqmwRESk2nu/7KdpAzfG9gk2O8ppVFgiIgLA3qPH+TkxkzsuC6FJA3ez45xGhSUiIgD8a+V+Gnm4cmf/dmZHOSMVloiIkJJTxA/bjjCub1uae3mYHeeMVFgiIsIHvx7AzdWFey4PMTvKWamwREScXMaxEr7ZnMYtMUH4Nm1gdpyzUmGJiDi5Gb8lYxhwnx3cQuRcVFgiIk4s+3gZ65JzufvyEAKbNzr/C0ykwhIRcWIzfjvA3qPHGdPb/q67+jMVloiIk8o5UcacdYcZ2S2AkFZeZsc5LxWWiIiT+mhVMmWVVfxtcJjZUWpEhSUi4oTyisr5z9pDXNe1De19Gpsdp0ZUWCIiTuijVcmUVFTxkIPsXYEKS0TE6eQXlfPp7ylc08WfMN8mZsepMRWWiIiTmbn6IMUVVTwcG252lAuiwhIRcSIFxeX8+/cUhnf2p4Of4+xdgQpLRMSpzFp9kBNllTwU6zjnrv6gwhIRcRLHSir4ZE0Kw6Ja07F1U7PjXDAVloiIk/hkzUGOl1U63LmrP6iwREScQGFpBbNWH2RoJz86tXG8vStQYYmIOIV5G1Px8nRz2L0rUGGJiNR7x4oreGt5EtEB3nQO8DY7zkVTYYmI1HMfr07meGklj8Z1MDvKJVFhidO466678PX1pXPnztXL8vLyiIuLIzw8nLi4OPLz86sfmzx5MmFhYURERLB06VIzIotcsryicmatPsg1XfyJ9HfMc1d/UGGJ07jjjjtYsmTJKcumTJlCbGwsSUlJxMbGMmXKFAASExOZO3cuu3btYsmSJTz44INUVVWZEVvkknz46wGKK6p4dIjjnrv6gwpLnMbAgQNp0aLFKcsWLFhAfHw8APHx8cyfP796+ZgxY/D09CQkJISwsDA2bNhg68gilyTreCmz16Ywsmsbwh1sVoszUWGJU8vMzMTf3x8Af39/srKyAEhPTycoKKj6eYGBgaSnp5uSUeRiTV95gIoqg0eGOPa5qz+4mR1AxB4ZhnHaMovFcsbnzpgxgxkzZgCQnZ1dp7lEairjWAmfrT/MqB6OcTfhmtAeljg1Pz8/MjIyAMjIyMDX1xc4uUeVmppa/by0tDTatGlzxq8xfvx4EhISSEhIwMfHp+5Di9TA+7/sxzAMHhrs+Oeu/qDCEqc2YsQIZs+eDcDs2bMZOXJk9fK5c+dSVlbGwYMHSUpKonfv3mZGFamxtPxivtyYyuiYIIJaNDI7Tq3RIUFxGrfeeisrV64kJyeHwMBAnn/+eZ5++mlGjx7NzJkzCQ4OZt68eQBERUUxevRoOnXqhJubG++//z6urq4m/wtEaubd5fuxWCz8zYHuJlwTljMdq/8f53xQRE4VExNDQkKC2THEiaXkFBH7xq/c1rctk0ZEmR2nJs58cvgMdEhQRKQeeWd5Eu6uFh4c1N7sKLVOhSUiUk/szzzOwh0ZxPdrh2+TBmbHqXUqLBGRemLaz/sIbNaQ8QNDzY5SJ1RYIiL1wLbUAhbvPMqIbm1o2djT7Dh1QoUlIlIPTF26hxZeHtwzoH7uXYEKS0TE4a1OymHN/lz+NiiMxp7192olFZaIiAMzDIOpS/cQ0KwhY/sGmx2nTqmwREQc2JKdR9medoxHh4Tj6Va/L25XYYmIOKjKKiuv/bSXcN/G3Ngj0Ow4dU6FJSLioL7dnE5ydhFPXBWBq0uNJ4xwWCosEREHVFpRxZvL9tEtqBlDO/mZHccmVFgiIg5ozrpDZBwr5clhEWe9V1t9o8ISEXEwx0sreP+X/QwIb0X/9q3MjmMzKiwREQczZ91hWjX25MmrOpodxaZUWCIiDiSzsJR3lifRwa8xXQK9zY5jUyosEREH8tayfVRarTw1LNLsKDanwhIRcRD7Mo/z5cZUbuvbjuCWjcyOY3MqLBERB/Hq4j14ebrx0OAws6OYQoUlIuIA1h7IZfmeLB68MozmXh5mxzGFCktExM5ZrQaTF++mjXcD7rysndlxTKPCEhGxcz/uyGB72jEeHxpBA/f6PcHtuaiwRETsWFllFVOX7iHSvynXdw8wO46pVFgiInZszrrDpOaVMPHqjk4xwe25qLBEROzUsZIK3l2RxIDwVgzs4GN2HNOpsERE7NT0lfuJ8GvC08Ocawqms1FhiYjYoUO5RcxcfZDgFo2ICnCuKZjORoUlImKHJi/ag7urCxOuijA7it1QYYkAb775JlFRUXTu3Jlbb72V0tJS8vLyiIuLIzw8nLi4OPLz882OKU5iXXIuS3Yd5YEr2uPbtIHZceyGCkucXnp6Ou+88w4JCQns3LmTqqoq5s6dy5QpU4iNjSUpKYnY2FimTJlidlRxAlarwUs/JtLGuwH3Dgw1O45dUWGJAJWVlZSUlFBZWUlxcTFt2rRhwYIFxMfHAxAfH8/8+fPNDSlO4ZvNaexML+TJYR2d+iLhM1FhidMLCAjgiSeeIDg4GH9/f7y9vRk6dCiZmZn4+/sD4O/vT1ZW1hlfP2PGDGJiYoiJiSE7O9uW0aWeKSqr5LWle+ka1IwRXduYHcfuqLDE6eXn57NgwQIOHjzIkSNHKCoqYs6cOTV+/fjx40lISCAhIQEfH10rIxfvw18PkHW8jP93bSQuTn6R8JmosMTpLVu2jJCQEHx8fHB3d+fGG2/k999/x8/Pj4yMDAAyMjLw9fU1OanUZ0cKSpixKplro/3p2baF2XHskgpLnF5wcDDr1q2juLgYwzBYvnw5kZGRjBgxgtmzZwMwe/ZsRo4caXJSqc+mLtmD1YCnr9ZFwmfjZnYAEbP16dOHm266iR49euDm5kb37t0ZP348J06cYPTo0cycOZPg4GDmzZtndlSpp7YczmdH+jHuuTyEwObOdyfhmrIYhnGux8/5oIicKiYmhoSEBLNjiAOxWg1umP47RwtKWPbYFTRp6G52JFur8ck6HRIUETHRN5vT2JZawFNXd3TGsrogKiwREZMcL63g1SV76RHcjOu7Ofe9rmpC57BEREzy7or95BaVMeuOGA1jrwHtYYmImGB/1glmrT7ILTFBRAc2MzuOQ1BhiYjYmGEYvLgwkYYerjyh2dhrTIUlImJjK/Zk8eu+bB4d0oFWjT3NjuMwVFgiIjZUVlnFCwsTCfNtzO392podx6GosEREbOizdYexWg2eu64T7q56C74Q2loiIjaSll/M1KV76BzgzYBwTZR8oVRYIiI28sIPiViw8M9rO5kdxSGpsEREbOCXPVn8lJjJQ7FhBDRraHYch6TCEhGpY6UVVUz6YRehPl7cc7lue3+xNNOFiEgd++DXAxzKLWbO3X3wcNN+wsXSlhMRqUOHc4v518oDXBPtz+XhrcyO49BUWCIidcQwDCb9sAt3FwvPXqOBFpdKhSUiUkeW7c5ixZ4sHh3SgdbeDcyO4/BUWCIidaC4vJIvNhwmpm1z7risndlx6gUVlohIHXh7eRIr9mQx4aoIzWhRS7QVRURq2e6MQj5edZAxvYLoE9rS7Dj1hgpLRKQWWa0GE7/dQbOG7jx9dUez49QrNimszYfzeW3pHlusSkTEVJ9tOMzW1AKevbYTzRp5mB2nXrFJYa1OyuH9Xw6wYk+mLVYnImKKrMJSpi7ew+VhrRjZrY3ZceodmxTW/Ve0J9y3Mc/O30VRWaUtVikiYnPPL0ykrMrKi9d3xmKxmB2n3rFJYXm4uTBlVBeOHCth2k/7bLFKERGb+mVPFj9uz+ChQWGEtPIyO069ZLNBFz3btmBcn7b8+/eDbEstsNVqRUTqXHFZJc8u2EmYb2PGX6HJbeuKTUcJThgWgU8TT57+dgcVVVZbrlpEpM688fNe/Jp48vL1nfF0czU7Tr1l08Jq2sCd50d0ZndGITNXH7TlqkXOqaCggJtuuomOHTsSGRnJ2rVrycvLIy4ujvDwcOLi4sjPzzc7ptihrakFzFqTQkf/prrmqo7Z/DqsYZ1bc1WUH28t28eh3CJbr17kjB555BGGDRvGnj172LZtG5GRkUyZMoXY2FiSkpKIjY1lypQpZscUO1NeaeWpr7fj26SBrrmyAVMuHH5+RGfcXVx4Z3kShmGYEUGkWmFhIb/99ht33303AB4eHjRr1owFCxYQHx8PQHx8PPPnzzcxpdij6SsPsDfzOC9d35kmDdzNjlPvmVJYrb0bMGlkFN9sTuerhFQzIohUS05OxsfHhzvvvJPu3btzzz33UFRURGZmJv7+/gD4+/uTlZV1xtfPmDGDmJgYYmJiyM7OtmV0MdG+zOO890sSI7q2YUgnP7PjOAXTpma6oVsAfUNb8NLC3Rw9VmpWDBEqKyvZvHkzDzzwAFu2bMHLy+uCDv+NHz+ehIQEEhIS8PHxqcOkYi+qrAZPfr2dxp5uPHed7nNlK6YVlouLhVdHRVNhtfKP73bo0KCYJjAwkMDAQPr06QPATTfdxObNm/Hz8yMjIwOAjIwMfH19zYwpduTfv6ewNbWASSOiaNnY0+w4TsPUyW/btvRiwlUdWb4ni/lb082MIk6sdevWBAUFsXfvXgCWL19Op06dGDFiBLNnzwZg9uzZjBw50syYYicO5xbz+tK9DIrwYURXTb9kS25mB7ijfzsW7chg0veJXBbWCt8muiun2N67777L2LFjKS8vJzQ0lE8++QSr1cro0aOZOXMmwcHBzJs3z+yYYjLDMJj43XZcLPDyDV00/ZKNWc5zKM4mx+n2Z51g+DurGBzhy/RxPfRDIA4rJiaGhIQEs2NIHflq42HmrD/MzTFB3Na3rdlx6osav+Hbxf2wwnwb8/chHViy6yiLdhw1O46IyGnS8ot5YeFuGrq7MLZ3kNlxnJJdFBbAvQNCiA705v8t2EnuiTKz44iIVLP+d1SgYRi8fnM3XFzs5q3TqdjNVndzdWHqTdFUWa3865cDZscREak2Z/0hfj+Qyz+u6URQi0Zmx3FadlNYAB1bN+X+K9ozc81BFm4/YnYcERFScoqYvGgPAzv4cKsOBZrKrgoL4J4BoXQNasY/5+8kq1AXFIuIeaqsBhO+3oabq4VXR2lUoNnsrrDcXF2YdnNXSsqrePpbXVAsIub5ZM1BNqbk89x1Ufh7NzQ7jtOzu8KCk6MGnxrWkRV7sjTXoIiYYn/WCaYu3cuQSF9G9QgwO45gp4UFJy8o7hfakhd+SCQ1r9jsOCLiRCqrrDw+bxuNPFx55UYdCrQXdltYLi4WXrs5GovFwhPztmG16tCgiNjGp2sPsS21gBdGdtbsO3bEbgsLILB5I/7ftZ1YfzCPT35PMTuOiDiBnenHeGVRInddFsJ10f5mx5H/YdeFBXBzTCCxHX35cuNhko4eNzuOiNRjJeVVPDx3C60aN+Dh2DAdCrQzdl9YFouFKTdGc6ykgke+3EpZZZXZkUSknnp5USLJ2UVMG92VZo08zI4jf2L3hQXg09STl6/vQmJGIdN+2md2HBGph1bsyWTOusPcOyCEy8JamR1HzsAhCgtgSCc/xvYJZsZvyazZn2N2HBGpR3JOlPHk19vp2LoJT1wVYXYcOQuHKSyAf17TiVAfLx7/ahv5ReVmxxGResAwTk5sW1hayTu3dsfTzdXsSHIWDlVYDT1ceWdMd3KLynjmO82CISKX7rP1h1mxJ4uJV3ekg18Ts+PIOThUYQF0DvDmiaERLN55lHkJaWbHEREHtj/rBC/9mMiA8FbE92tndhw5D4crLIB7B4TSv31LJv2wi4M5RWbHEREHVFZZxTvLk2jk4ca0m7vi4qIh7PbOIQvLxcXCtNFd8XC18PayfZRXWs2OJCIOZsriPXy/7QjTbo7Gt6lms3AEDllYAP7eDXn95m7M33qEqUv2mB1HRBzIz4mZfLImhTv6t2NQRz+z40gNOWxhwcmh7rf3a8vHqw+yLDHT7Dgi4gDSC0p4Yt42oto0ZeLwjmbHkQvg0IUF8MzwSDr5N+WJr7dxpKDE7DgiYscqq6w88sUWKqusvPeXHhrC7mAcvrAauLvy3l+6U1Fp5ZG5J38QRUTO5M1l+0g4lM8rN3YhpJWX2XHkAjl8YQGE+jTm5Ru6sDEln7eWJZkdR0Ts0KqkbP618gCjYwIZ2U03ZHRE9aKwAK7vHsDomEDeX7mf1UmauklE/k/W8VL+/uVWwnwaM2lElNlx5CLVm8ICmDQiijCfxjz65VYyC3U+S0TAajV47MttHC+t5L2/9KCRh5vZkeQi1avCauThxvtjexDm68WTX+/Q+Sy5IFVVVXTv3p1rr70WgLy8POLi4ggPDycuLo78/HyTE8rFmL5yPwUl5UwaEUVEa0295MjqVWEBdPBrwk09g/h1XzbTftatSKTm3n77bSIjI6v/PmXKFGJjY0lKSiI2NpYpU6aYmE4uxi97s3j9532E+zZmTK8gs+PIJap3hQVwU89Abu0dzPSVB/hp11Gz44gDSEtL48cff+See+6pXrZgwQLi4+MBiI+PZ/78+Salk4uRmlfMo3O3EuHXhFduiNbdg+uBellYAM9d14kuAd48Pm8bh3I136Cc26OPPsrUqVNxcfm/X4nMzEz8/f0B8Pf3Jysry6x4coFKK6p44LNNWA2DD8b1pKGHrreqD+ptYTVwd+VfY3vgYrFw/5zNlFZUmR1J7NTChQvx9fWlZ8+eF/X6GTNmEBMTQ0xMDNnZ2bWcTi6UYRg8O38nO9MLeXN0N9rpeqt6o94WFkBQi0a8NaYbe44W8s/5O3X/LDmjNWvW8P3339OuXTvGjBnDihUrGDduHH5+fmRkZACQkZGBr6/vGV8/fvx4EhISSEhIwMfHx5bR5Qzmbkxl3qY0HhocxpBOmiewPqnXhQUwKMKXhwaH8/WmNL7cmGp2HLFDkydPJi0tjZSUFObOncvgwYOZM2cOI0aMYPbs2QDMnj2bkSNHmpxUzmdbagHPLdjFgPBWPDqkg9lxpJY5xQUJj8SGc/RYCdN+3keH1k3oEdzc7EjiAJ5++mlGjx7NzJkzCQ4OZt68eWZHknPIKyrngTmb8GniyTtjuuOq+1vVO5bzHCarN8fQCorLGfHeGkorqvjhocvx0/1vpA7ExMSQkJBgdgynU2U1iJ+1gQ0peXx9fz+iA5uZHUlqrsafLOr9IcE/NGvkwUe3x3CirJL7/rOJskoNwhCpL6b9tJfV+3N4cWSUyqoec5rCAoho3YQ3Rndla2oBz2oQhki9sGhHBr/vz+X2fm25pVew2XGkDjlVYQEM6+zPw4PD+CohjU/XHjI7johcgp3px3jsq624uMA/rok8/wvEoTldYQE8OqQDQyJ9eWFhImsP5JodR0QuQlZhKfd+mkCLRh58eFuMbsboBJyysFxcLLx5SzfatWzEXz/fTFp+sdmRROQClFZUMf4/mygoruCj+Bh8mniaHUlswCkLC6BJA3c+uj0Gv6aePPXNdk6UVZodSURqwDAMJn67g62pBbx5S1ei2nibHUlsxGkLC07eqfjpYR1Zl5zHw19socqqQRgi9u6DX5P5bks6TwztwLDO/mbHERty6sICuCLCl0kjolixJ4uXfkw0O46InMPPiZlMXbqH67q24a+DwsyOIzbmFDNdnM9tfduSklPEzNUHCW3lxW392pkdSUT+ZHdGIY/M3UKXAG9eu0m3C3FGKqz/emZ4JIdyi5j0QyJBLRpxZcSZJzoVEdvLOl7KtKV78W548txzA3eNCHRGTn9I8A+uLhbeHtOdCL8m/O3zLew9etzsSCICFJVVcte/N7LmQC4f3d5T06o5MRXW//DydGPmHTE08nDlrn9vJPt4mdmRRJxaZZWVv32+mcQjhbw/tjudA5qZHUlMpML6E3/vhsyM70XuiTJeWbSb4nINdxcxg2EYPLtgF7/szebF6zszuKPubeXsVFhn0CXQmw9u68mCrek89PkWKqusZkcScTrTfz3AFxsO88CV7Rnbp63ZccQOqLDO4soIX14Y2Znle7J0t2IRG1uwNZ2pS/YyomsbJgyNMDuO2AmNEjyHcX3bkllYyrsr9uPXtAF/j9MdTEXq2toDuTwxbxt9Qlrw2s3RuOhGjPJfKqzzeCyuA0ePlfL28iRaezfg1t66fYFIXUk8Usj0lftp19KLGZrQVv5EhXUeFouFV27sQvaJMv7x3Q5aNfYgrlNrs2OJ1DuHcou4fdYG3F0tfP1Af7wbuZsdSeyMzmHVgLurC+//pQeDInx5aeFu1ifrliQitSmrsJRxM9dTZbXyn7t7E9CsodmRxA6psGrIy9ON127uirubC3fPTmBH2jGzI4nUC8eKK7ht5gZyT5TzyZ29CfNtYnYksVMqrAvQwsuDOXf3oVkjd+I/2cD+LM2GIXIpissruWv2Rg7mFDHjthi6BTUzO5LYMRXWBWrt3YA5d/fBxWJh3McbSM3TzR9FLkZ5pZUH5mxmy+F83h7TjcvDW5kdSeycCusitGvlxZx7elNSUcW4mevJKiw1O5KIQ6mssvLiwkQ2HcrnlRu6cHUX3ddKzk+FdZE6tm7KJ3f2Ivt4GU99s538onKzI4k4hCqrwYSvt/OfdYd4clgHxuhSEakhFdYl6BHcnI/jY9iWdoyxH6+noFilJXIuVqvBxG+3892WdCZcFcHt/ULMjiQORIV1ifq3b8Ubo7uyP+sE42au51hxhdmRROzSyclsd/JVQhoPDw7THYPlgqmwasGVEb58eFtP9h09wW2z1nOsRKXlSFJTUxk0aBCRkZFERUXx9ttvA5CXl0dcXBzh4eHExcWRn59vclLHZRgGz/+QyGfrD3PfFaGa5kwuigqrlgzq6Mv0cT3YnVHI7bM2UFiq0nIUbm5uTJs2jd27d7Nu3Tref/99EhMTmTJlCrGxsSQlJREbG8uUKVPMjuqQDMPg/V/28+/fU7jzsnY8Payjbm8vF0WFVYtiI/14/y892JV+jPhZGziu0nII/v7+9OjRA4AmTZoQGRlJeno6CxYsID4+HoD4+Hjmz59vYkrHZBgGLyxM5PWf9vHokHD+37WdVFZy0VRYtWxoVGve+0sPdqQd47kFu7Sn5WBSUlLYsmULffr0ITMzE3//k8Ot/f39ycrKMjmdY7FaDf45fyefrEnh7stDeCQ2XGUll0SFVQeGdW7NjNt68v22I4z9aL2GvDuIEydOMGrUKN566y2aNm1a49fNmDGDmJgYYmJiyM7OrsOEjqPKavD0t9v5bP3JGzD+85pIlZVcMhVWHRkc6ceM23uyN/M4Y2asI+u4Li62ZxUVFYwaNYqxY8dy4403AuDn50dGRgYAGRkZ+Pr6nvG148ePJyEhgYSEBHx8fGyW2V5VVll5Yt42vkpI45HYcJ68KkJlJbVChVWHBnf045M7enE4r5hbPlzHkYISsyPJGRiGwd13301kZCSPPfZY9fIRI0Ywe/ZsAGbPns3IkSPNiugwyiuqeOTLrdXXWf09roPKSmqN5Ty3ftd94WtBQkoed36ykaYN3fn83j60belldiT5H6tXr2bAgAF06dIFF5eTn+FeeeUV+vTpw+jRozl8+DDBwcHMmzePFi1anPNrxcTEkJCQYIvYdqekvIoHP9tEUXkVQyJ9GT+wvdmRxDHU+BONCstGdqQd47aZ6+gW3Jynr+5Ix9Y1P0cijsNZC+tYSQX3zN5IwqF8XhgRxW392pkdSRxHjQtLhwRtpEugN/Pu78/ujEJGf7CWjSl5ZkcSqRVZx0sZM2MdW1MLePfW7iorqTMqLBsK92vCNw/0p1UTT8Z9vJ6fEzPNjiRySQ7nFnPzB2tJySliZnwvro1uY3YkqcdUWDYW2LwRX9/fn47+TbnvPwl8tTHV7EgiF2VH2jFun7Wepg3c+ezePgzsoBGSUrdUWCZo4eXB5/f04fJwH/61cj/vLE/iPOcSRezKij2ZjP5wLRVVBm+MjqZHcHOzI4kTUGGZxMvTjY9vj2FABx/e+HkfT8zbTnml1exYIuf1+frD3DM7gfa+Xnz31/6E+2kAkdiGm9kBnJmHmwsvjIiilZcnby7bx5GCEj4Y1xPvRu5mRxM5jdVq8PpPe/nXygMMivDhvb/0wMtTbyFiO9rDMpnFYuGRIeG8eUtXEg7lceP0NaTmFZsdS+QUJeVVPP3tdr7YcJhbewfz0e0xKiuxORWWnbiheyD/ubsPOSfK+cf8HRr2LnYj41gJoz9cy7xNaTw6JJxXbuiMm6veOsT29FNnR/qGtuS7B/tzpKCUv3y0js/XHzY7kji5zYfzGfHeGpKzT/Dx7THE9w/RVEtiGhWWnQn1acw3D/Snf/tWPPPdDv45f4cGY4gpvt2cxpgZ62jo7sp3f72M2Eg/syOJk1Nh2SHvhu7MuqMX910Rypx1hxk3cz05J8rMjiVOorLKyjvLk3jsq230CG7Ggr9eRge/JmbHElFh2StXFwsTr47k7THd2JZawA3vr2FbaoHZsaSeyz5exm0zN/DhrwcYPzCU/9zdh+ZeHmbHEgFUWHZvZLeA6pkxbv5gLf9Zd0gXGUudSEjJ49p3V7H5cD7Pj+zMM8MjcdfgCrEj+ml0AF0CvXl1VDT9w1ry7PydPDJ3KyfKKs2OJfWEYRh8vCq5+nzV/L9exk09A82OJXIaFZaDaOHlwaz4Xky4KoKF248w4r3V7D163OxY4uCOFVfw/A+JvPTjbgZ39OX7hy4n0l8zV4h90v2wHNDvB3J4+IutBDZvyM0xgfyld7CGGtsJR7of1vrkXP7+5Vbyisp5Zngkt/Vrq58jMYPuh1Wf9W/fikWPXI6/dwP+8d1O7v10E7kaRSg1VFFlZdpPe7n1o3V4uLnw5X39uL1/O5WV2D3tYTkwq9Vg1pqDTF2yl6YN3Xnt5mgGRfiaHcup2fse1qHcIh6Zu5WtqQXc3DOQSSOiNMWSmK3Gn5RUWPXA7oxCHp27lb2Zx4nv15aJwyNp4O5qdiynZK+FZRgG3287wjPf7sDFxcLkG7voZotiL2pcWPpoVQ9E+jdlwd8uY+qSvcxac5DU/BIeGhxGd92jSICjx0r5x3c72HI4nz6hLXnx+s4ENGtodiyRC6Y9rHrmt31ZTF26l8Qjhdx1WQiPD42goYf2tmzFnvawDMPgq4RUXlq4mwqrlQlXdeSO/u1wddG5KrEr2sNyVgM7+NI9uDlTFu/h49UH+SkxkymjutC/fSuzo4kNpeUXM/HbHaxKyqFPSAteHRVNu1ZeZscSuSTaw6rH1iXn8vQ320nJLeYvvYN5clgEzRppmp26ZPYeVmWVlU/XpvDp2kNkHS9j4tUdGdunLS7aqxL7pUEXclJJeRVvLdvH9rQC9maeYMJVEYyOCdJhoTpiZmFtOpTPP+fvZHdGIWN6BfHXQWEEtWhkShaRC6DCklMlHilk0ve72JCSR5cAbyaNiKJnWw3KqG1mFFbeiTKmLt3L3I2ptG7agOeu68Swzq11XZU4ChWWnO6Poc2TF+3haGEpN/YI4OlhHfFt2sDsaPWGLQurrLKK/6w9xJx1h8grLufWXsE8HBuu66rE0WimCzmdxWJhZLcAlj9+BQ9e2Z6F2zIYPO1XZq46SFllldnx7NKSJUuIiIggLCyMKVOmmB0HOPnBY9GODOLe+I2XftxN25ZefPtAfyYOj1RZSb2mPSwndjCniBcXJnK8tIL0/BIeGRLOqB6BuOmWEgBUVVXRoUMHfv75ZwIDA+nVqxdffPEFnTp1Outr6noPa/OhfF5ZtJuEQ/lE+DXhmWsiuaKDT52tT8QGNKxdzi+klRez7ujF6qRsXvtpH099s4MPf03m0bgOXNvF3+lHlm3YsIGwsDBCQ0MBGDNmDAsWLDhnYdWVnenHeGvZPjKOlZJZWMbkG7twc099uBDnosISLg/34bKwVizbncXrS/fy8Bdb+Ncv+3liaASxkb5Oe/I+PT2doKCg6r8HBgayfv16m2ZIPFLIW8v28VNiJt4N3bl3QAh39A+hcQP96orzOechwWHDhhk5OTm1sqLs7Gx8fBzv0IWj5oaLz15QXEFmYSlVVgMPNxdaNfbEu5F7zffba4E9bPf8/HwKCwtp27YtALm5uRQVFREcHHzK87Kzs/nj96SsrIxu3bpd8rpLyqvIKSqjoLgCF4sFn8aetGzigWsdfniwh21+MRw1Nzhu9trMvWnTpqWGYQyryXNtdg7L7AsqL5aj5oZLy15RZWXxzgzeWpZEcnYRAc0actflIdzSK4jGNjixbw/bfe3atUyaNImlS5cCMHnyZAAmTpx41td4eXlRVFR0UeuzWg1+2ZvFR6uSWZecxxUdWtE1qDl3Xx6Cd0P3i/qaF8IetvnFcNTc4LjZazm3zmHJpXF3dWFE1wCu7dKG5Xuy+Oi3ZF5cmMjby/Yxtm9b7uzfrt4Ph+/VqxdJSUkcPHiQgIAA5s6dy+eff17r6yktr2T+1iN8tCqZA9lFtPFuwD+GR3JL7yCaNqj7ohJxFCosOScXFwtxnfyI6+THlsP5fLQqmQ9/PcDPiZlEtG7CzT0DGRDuUy9nznBzc+O9997jqquuoqqqirvuuouoqKha+/r7s47zxYZUftuXzYHsE0T6N+XtMd0Y3sUfdw2mEDmNzQpr/PjxtlpVrXLU3FD72bsHN+dfY3tyKLeIH7YdYebqg/y4PYPWTRswqmcAN/UMIqSWJli1l+0+fPhwhg8fXuPnt2p17kmGS8qr+HFHBnM3HCbhUD7urhaGdmrNKzd0IaZdc1MHuNjLNr9QjpobHDe7Wbl1HZZctLLKKpbvzmJeQiq/7svGakCvds25uWcQw6P9bXKuy96c6dh+ZZWV3w/k8uP2DDam5JGcU0RoKy/G9A7ixh6BtGrsaVJaEbugqZnEtjILS/lmcxpfJ6SRnFNEB7/GtPdpTGykH4MifGjpJG/KfxRWZZWVDSl5LNyewZKdR8krKqexpxu39WvLFR186BPSwmkvFxD5E3OmZpo3bx5RUVG4uLic9ilz8uTJhIWFERERUT3q6s/y8vKIi4sjPDycuLg48vPzazNejdxyyy1069aNbt260a5du7MOUW7Xrh1dunShW7duxMTE2DbkWUyaNImAgIDq/IsWLTrj8+piuiG/pg148Mowlj9+Bd880I9ruviz6VA+T8zbRszLy7hp+u9MX3mApMzj/PlD0oQJE+jYsSPR0dHccMMNFBQUnHEd9rTNz7QNMwtLyS8u56EvtjDo9ZWMnbGWz9ckUZi0kSf7Nyfhn0N4alhH+oa2NKWsUlNTGTRoEJGRkURFRfH222+f9pyVK1fi7e1d/TP0wgsv2Dzn2Zzv+28YBg8//DBhYWFER0ezefNmE1Kebu/evdXbs1u3bjRt2pS33nrrlOfYy3a/66678PX1pXPnztXLavq+bJNpzAzDONd/FyQxMdHYs2ePccUVVxgbN26sXr5r1y4jOjraKC0tNZKTk43Q0FCjsrLytNdPmDDBmDx5smEYhjF58mTjySefvNAIteqxxx4znn/++TM+1rZtWyM7O9vGic7tueeeM1577bVzPqeystIIDQ01Dhw4YJSVlRnR0dHGrl276iRPVZXV2Jaab0z7aa8x/O3fjLZPLTTaPrXQGPDqCmPqkj3Gou1HjOzjpcbSpUuNiooKwzAM48knnzzr991etvkf23DX3n3Gil1HjPCbnzIGTl5qtH1qoeHROsyIeelnY/S0hUb/0Q8YJ0rLjbVr1xq9e/c2O7Zx5MgRY9OmTYZhGEZhYaERHh5+2vf+l19+Ma655hoz4p3X+b7/P/74ozFs2DDDarXazTb/s8rKSsPPz89ISUk5Zbm9bPdff/3V2LRpkxEVFVW9rCbvy5f4vnK+Hqr+r1ZPMkRGRp5x+YIFCxgzZgyenp6EhIQQFhbGhg0b6Nev32nPW7lyJQDx8fFceeWVvPrqq7UZscYMw+Crr75ixYoVpqy/rthyuiEXFwvRgc2IDmzGY3EdyDhWwvLdWazYncm8hFTeP14GQHsfLxIq99G2pRcBXfqxdvHXtZ7lUhmGweG8YnZnFPL92l24D5/IQz+kcyC7CJfQy6g8nsfEq3vx7uLGbHgmlvvvv5/rr78SL093+vbtS0FBARkZGfj7+5v2b/D3969ef5MmTYiMjCQ9Pd2UqabqwoIFC7j99tuxWCx2s83/bPny5bRv3776YnR7M3DgQFJSUk5ZVpP3ZVu9r9jkrHh6ejp9+/at/ntgYCDp6emnPS8zM7P6h8vf35+srCxbxDujVatW4efnR3h4+Bkft1gsDB06FIvFwn333Wc3o33ee+89Pv30U2JiYpg2bRrNm596zyszpxvy927IuL5tGde3LeWVVnYeOcaGg3msT87lp12ZpOQWA+6ERNzENe+s4rKwVjRv5EG4X2OCmjcEN4863+aVlVYyCks5lFtMUtZxkrJOcCi3iO1pxzheWklIKy9Sc0tp7OHOkEg/nr22JQfWLWXLxvXcd8UYPnJ3xWKxnHE7p6en282bZ0pKClu2bKFPnz6nPbZ27Vq6du1KmzZteP3112t1KP+lON/vnL1vc4C5c+dy6623nvExe93uNXlfttX7ygUX1pAhQzh69Ohpy19++WVGjhx5xtcYZxjYYeYJ55r8G7744ouz/mABrFmzhjZt2pCVlUVcXBwdO3Zk4MCBdZb5D+fK/sADD/Dss89isVh49tlnefzxx5k1a9YpzzPre3Gu3J/cORLDMHj6pddZn3SUUXfcz4aUfLalFrD+YB4tvTxwc7HATW8SHeFDcUkZU9euYxs/E9YuGJ8mnjRp6E4jD1c83Vxo4OaGq8vJi58NDKwGVFRaKa+yUlpppbS8iqLySgpLKkjLL8HVxcK2tGN4urmw52gheUUVtGvRiJS8YrwbutPBrzGjY4II921MpzZN2bXmZ5b/vIqJw+MBSN1w+ja0t5/5/3XixAlGjRrFW2+9RdOmTU95rEePHhw6dIjGjRuzaNEirr/+epKSkkxKeqrz/c7Z8zYHKC8v5/vvv6+eMeV/2fN2rwlbbfsLLqxly5Zd8EoCAwNJTU2t/ntaWhpt2rQ57Xl+fn7Vu/AZGRn4+vpe8Lpq4nz/hsrKSr799ls2bdp01uf8kd/X15cbbriBDRs22KSwarr97733Xq699trTltf0e1Hbzpf7008/5bdF37J8+XIaNfq/27oXllaQnHWCIwWlJOecoLTCSnpBCUmt27Eh5RirjqTQtmUjGnq4UVpRhburBTcXF8qrrLhaLLi6gIEFq2FgtRrVFzhbgEqrwdFjpUQHNcPd1ULnAG+6BXnTtqUXbVt6EerjhU9jz9N+8YqCz78NzdrO51NRUcGoUaMYO3YsN95442mP/2+BDR8+nAcffJCcnJzzXl9mC+f7nbPXbf6HxYsX06NHD/z8/E57zJ63e03el2227c9zkuui/HnQxc6dO08ZdBESEnLGQRdPPPHEKSf3JkyYcLERLsnixYuNgQMHnvXxEydOGIWFhdV/7tevn7F48WJbxTurI0eOVP/5jTfeMG655ZbTnlNRUWGEhIQYycnJ1SdHd+7cacuYp1m8eLERGRlpZGVlnfU5Z9vmlZVVRv6JMuNwbpGReKTA2JySZ6xLzjFWJWUZv+3LMn7be/L/v+/PNhJSco2daQXG/qxCI6uw1CivOP1nsCbOtQ179uxpGIZhLFy48JQBAL169bqoddUmq9Vq3HbbbcYjjzxy1udkZGQYVqvVMAzDWL9+vREUFFT9dzPV5HfOHrf5/7rllluMWbNmnfExe9ruBw8ePGXQRU3ely/xfaXGgy5qtbC+/fZbIyAgwPDw8DB8fX2NoUOHVj/20ksvGaGhoUaHDh2MRYsWVS+/++67q8stJyfHGDx4sBEWFmYMHjzYyM3NvdAItSI+Pt6YPn36KcvS09ONq6++2jAMwzhw4IARHR1tREdHG506dTJeeuklM2KeZty4cUbnzp2NLl26GNddd111gf1vdsM4OZoqPDzcCA0NtYvs7du3NwIDA42uXbsaXbt2Ne677z7DMOx7m59pG06fPt0IDg42DONkOTz44INGaGio0blz51M+wJll1apVBmB06dKlelv/+OOPxvTp06t/3t99912jU6dORnR0tNGnTx9jzZo1Jqc+6Wzf///Nbo/b/A9FRUVGixYtjIKCgupl9rjdx4wZY7Ru3dpwc3MzAgICjI8//vis78u1+L5S48LShcMitchRZ98WMZE5Fw6LiIjUFRWWiIg4BBWWiIg4BBWWiIg4BBWWiIg4BBWWiIg4BBWWiIg4BBWWiIjYzNSpU3nnnXcA+Pvf/47FYlkBYLFYYi0Wy5xzvfZ8Fw6LyAWwWCxLDMMYZnYOEXtlsVj6Ao8bhnGzxWJZBXgClwHPAEcNw/jwbK/VHpZILVJZiZzXJqCnxWJpApQBa4EYYACw6lwvtMn9sERERAAMw6iwWCwpwJ3A78B2YBDQHth9rtdqD0tERGztN+CJ//5/FXA/sNU4zzkqFZaIiNjaKsAfWGsYRiZQynkOB4IGXYiIiIPQHpaIiDgEFZaIiDgEFZaIiDgEFZaIiDgEFZaIiDgEFZaIiDgEFZaIiDgEFZaIiDiE/w+nNHQACEDmawAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sympy.plotting.plot(f);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A neat parabola. We can see from the plot that the minimum of $f(w)$ is reached somewhere between $w=-2.5$ and $w=0$. SymPy can tell us the derivative, and the value where it is zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADwAAAAQCAYAAABKkhw/AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACn0lEQVRIDdWX61EbQQyAD8YF8OgAOgCnA9MBhA6ggzD8sv9loAOSChjoIO4gQAfQQTzuwPk+ze6xdzmG4IwPohkhraRdPVbaw2uLxaKaTCY7VVWdgcIQnIFnyB8U/E+QcjlNMW9Aze0C+VTZ2ng8VnCF4ECBAH8B+QIeZEPlHx2I1QRNLidcwR8iuwGP4G/XNQBrA3iNvO05qGGvgO8ReLKkU/edsN8kM8TNsjhXYMIj8BEjq1OChhvI7YA+wTjasfytf0fQixIDiL/mFQxAE9trK1QmWNZ53t8bJQdz2SwdIsu3faV8gOCoNCj4PXn08XBBTdyW2QbtiG/QAHhtT6ExGlCdOCr78I0Kx4ae/uDb7o2RhY94bek/AKUJ2Mr55dbmHPkl9CfoISU4H2Xrm6TrYWnUF2/8oI+uF+CF3WXfnQmj9LG6TQlWHsDaRAVf81lwz3+s5I+8xN7W0tFTlvVJ8f8AXoJ27zV4Dx+t7Qw3AIW9/pSMs851tDaCz+DXrEBuMWx3kyzhDt2LCaPTj4Vqw5YC9I0vRzIykZdGsH1OrLH34uy4G+hmI2EEzugWtP4mu4t1zCHUKplcPb/wBj1HlwvCMuAxM10U+66EquRjB+r4vAnYY/E9ox2LLW2co7qlMTKZXWhdQXgdl7NpIaxyFABeUNa4XfQ6bjvVdtVwjwPb10vphEg4BfgJWj5SbrAIs2KnybfbtEt2zFmNIhRnrJL1InILl36GaTEdEJgB+0hN4eNbVViOkJWt5a1F22iDzjZx5jwjIMnqByzLe6LtCzPGPIZ+Nuf+L+2s1QG3ArN990sZ61wUq/kLdJ6/g3aCssY3mvWbIAW41AzriP1eQj2W8Ob2/OPBX0sfCf414ddyWX/N4B30dom4EvgNCwYAZwDIeMAAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle 2 w + 3$"
      ],
      "text/plain": [
       "2⋅w + 3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fprime = f.diff(w)\n",
    "fprime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACsAAAAzCAYAAAAO2PE2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACw0lEQVRoBe2a31HjMBCHk5srIFwJoQP+dBA6IC1ACcw9Ja/QAbQQOgA6IHQAJcB1cPf9fBYjyZYszTBrHrwzGknrtfV5vVopmsw3m81i1iPb7fZPj9pExdi9TD8YfU/5iMq1CVV6kF3EI77dHM++0rjkbR7T9457BbYLCNY/LTEYdMl4l+2Y+tTqX5c6ygwWIMFdUTvYGe1zdA/Ua8o97awoZq1En/KiBXRjutD77RS52hL2BRBlGJVGAP9sO12utgwDefHAh/G8fOvrU21LzwYMgK5QKEUqE90FFxMdM8+68QE7oi3QU4pC45lSJGPAClBlBriywZ7622UD8QUCpNKVJtmO9iK42NMxi1lgjlR6GFwYKDSyYgYLhfYg+uSDHkwRW8asPvcjsKp9OWk7boHwrwVtS9irYGQ6gGuCydNKX/FLxOYzM1hg7igrir8ALCE6QzfoVZGbwWqwFqoITPaxWE6weOzq/gRb7bLCGybPFjqq2mzybLXLCm/I5lnyolaXJ4rqUtF2r9kCuhvo/3XtoRrbecpmCFZL4HHq5lJ9DqD0GbKbYrbGWzW2k2drvFVjm51gNQ8qsWWiLbFz+1ptut/VRx9kj9SzzGBb0FvqMwdDW+cG+qlTtKe1jNnmQMOBqgZSXlZ61HnsoFjCrqB5BTBeYLQZX6BXiGTFElZQb0DJk30Sv0THxjJm153R/yuaswReYnCSWXq2wwqgQPX5XYbo2PiKUWEB0cS6B/rGh0q1R4MFUD/JFcOp8OgwjwILoI7sf1F/5twOWY/CHBZAncIc+h6lvVTp4QtUprAAaUKdUscTSi+gpTcrlqlLntOE0uGcf4QkQB0rDU4yM1iAHigCVrzGMphjdYMZLJ47jAlr+6YxWwsX20+wsUe+qu9iVnku2PXQT+2Ovmrs5HNiFgwbNhcGSiX+HyS+258jGp5/R7zM+TmH0KYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left[ - \\frac{3}{2}\\right]$"
      ],
      "text/plain": [
       "[-3/2]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.solve(fprime, w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks about right: $-3/2$ or $-1.5$. \n",
    "We could have also solved this by hand, because it's a simple function. \n",
    "But for more complicated functions, finding the minimum analytically could be more difficult. \n",
    "Instead, we can use the iterative method of gradient descent. \n",
    "\n",
    "The idea in gradient descent is to find the value of $w$ at the function minimum by starting with an initial guess, then iteratively taking small steps down the slope of the function, i.e., in the negative gradient direction. \n",
    "To illustrate the process, we turn the symbolic expression `fprime` into a Python function that we can call, and use it in a simple loop taking small steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "function"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fpnum = sympy.lambdify(w, fprime)\n",
    "type(fpnum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yep. We got a Python function with the [`sympy.lambdify()`](https://docs.sympy.org/latest/modules/utilities/lambdify.html) method, whose return value is of type `function`. Now, you can pick any starting guess, say $w=10$, and advance in a loop taking steps of size $0.01$ (a choice we make; more on this later):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.4999999806458753\n"
     ]
    }
   ],
   "source": [
    "w = 10.0   # starting guess for the min\n",
    "\n",
    "for i in range(1000):\n",
    "    w = w - fpnum(w)*0.01 # with 0.01 the step size\n",
    "\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That gave a result very close to the true value $-1.5$, and all we needed was a function for the derivative of $f(w)$. This is how you find the argument of the minimum of a function iteratively. \n",
    "\n",
    "##### Note\n",
    "\n",
    "> Implied in this method is that the function is differentiable, and that we can step *down* the slope, meaning its second derivative is positive, or the function is _convex_.\n",
    "\n",
    "<img src=\"../images/descent.png\" style=\"width: 400px;\"/> \n",
    "\n",
    "#### Gradient descent steps in the direction of the negative slope to approach the minimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression\n",
    "\n",
    "Suppose you have data consisting of one independent variable and one dependent variable, and when you plot the data it seems to noisily follow a trend line. \n",
    "To build a model with this data, you assume the relationship is _linear_, and seek to find the line's slope and $y$-intercept (the model parameters) that best fit the data. \n",
    "\n",
    "Though this sounds straightforward, some key ideas of machine learning are contained:\n",
    "\n",
    "- we don't _know_ the true relationship between the variables, we _assume_ it is linear (and go for it!)\n",
    "- the model we chose (linear) has some parameters (slope, intercept) that are unknown\n",
    "- we will need some data (observational, experimental) of the dependent and independent variables\n",
    "- we find the model parameters by fitting the \"best\" line to the data\n",
    "- the model with its parameters can then be used to make _predictions_\n",
    "\n",
    "Let's make some synthetic data to play with, following the example in Eric Ma's tutorial [1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAX70lEQVR4nO3db4xcV3nH8d/jzYoMgmYdxaTOEuOoQg6QtLZYISRLLXFIk0JKDLywqEojUdV9AVIchVQLvCBVVWE1lPCiVSu3oKYqalMpIYkSqtSN3SKsNmWdmIDlpEEF2mzcxJSsAHmVru2nL2bGzM7eM3Pu37l37vcjRbs7O7tzVpF+c/yc55xj7i4AQPNsmvQAAADZEOAA0FAEOAA0FAEOAA1FgANAQ11S5YtdccUVvn379ipfEgAa7/jx4z909y3Dj1ca4Nu3b9fS0lKVLwkAjWdmP0h6nBIKADQUAQ4ADUWAA0BDEeAA0FAEOAA0VKVdKADQNg8/s6x7n3heL62s6qq5ju6+eYf27pov5HcT4ABQkoefWdanHvq2VtfOS5KWV1b1qYe+LUmFhDglFAAoyb1PPH8xvPtW187r3ieeL+T3E+AAUJKXVlZTPZ4WAQ4AJblqrpPq8bQIcAAoyd0371BndmbdY53ZGd19845Cfj+LmABQkv5CJV0oANBAe3fNFxbYwyihAEBDEeAA0FAEOAA0FAEOAA1FgANAQxHgANBQBDgANBR94ACgco99LQsBDqD1yj72tSyUUAC0XtnHvpaFAAfQemUf+1oWAhxA65V97GtZxga4mV1qZv9uZt8ys5Nm9vu9xy83s8Nm9kLv4+byhwsAxSv72NeyxMzAX5O0x91/SdJOSbeY2bslLUp60t3fKunJ3tcA0Dh7d83rcx+6XvNzHZmkuc6sLp3dpDsfOKHdB4/o4WeWJz3ERGMD3Lt+2vtytvefS7pN0v29x++XtLeMAQJAFfbumtexxT26b99OvXbugl49uybXzzpS6hjiUTVwM5sxsxOSXpF02N2fknSlu5+WpN7HNwV+dr+ZLZnZ0pkzZwoaNgCUo0kdKVF94O5+XtJOM5uT9FUzuy72Bdz9kKRDkrSwsOBZBgkAVYntSKnDxp9UG3ncfcXM/lnSLZJeNrOt7n7azLaqOzsHgFobDN7LOrMyk1bOrl0M4avmOlpOCPHBjpS6bPyJ6ULZ0pt5y8w6kt4r6TlJj0q6vfe02yU9UtIYAaAQ/eBdXlmVS1pZXdtQ677h2i1jO1JCZZYDFS96xtTAt0o6ambPSvqmujXwxyQdlHSTmb0g6abe1wBQW0nBO2h17byOPndmXUfK/FxHn/vQ9etm1qM2+FS56Dm2hOLuz0ralfD4/0q6sYxBAUAZYnZWvrSyGryIuF9+GbeY11/0LLucwk5MAK0Rs7My9JzB8kuMKrbhE+AAWiNpx+WgUbsvx5VfhlWxDZ/jZAG0Rr+kMaoLJVT2GDWj7szOrAv3qrbhE+AAWiVU3x4n1F443wv+SfSEE+AAEOHum3es6/2WfjbTzvqmkBcBDgARhssvdbh2jQAHgEgxM+0qt9gT4ACmXlWhWvUWewIcQKOkDeMqQ3XUSYYEOIBWyxLGeUM1zRtG1XdrspEHQGNkOas7T6gOH3417pyTqu/WJMABNEaWMM4TqmnfMKq+W5MAB9AYWcI4T6imfcMYvlsz6STDIlEDB9AYozbThOTp34653CHp9arqDSfAATRG1jDOGqpZ3jCqRIADKE1MB0fatsAqZ7h13H05iAAHUIqYlr+ie7TL2LAzqXNOYhDgAHJLCs6Y/usiN77U5aLhKtGFAiCXUK906OaawQ6OIje+ZOkRbzoCHEAuoeCcMUt8/mAHR5EbX6reBVkHBDiAXEIBed59bP91kRtfqt4FWQcEOIBcRgXk6y7ZpM2vnw1uaili48vDzyxr98EjWl5Z1fCc39Qt6ew+eCS4/b3JzN0re7GFhQVfWlqq7PUAlG948XBYZ3amtN2ISa9tknzgYxXjKJuZHXf3heHHmYEDyGVwFp2kzIXEpPq7S5ox0/DUdBoXNGkjBJBbv1f6msXHNwSntLFOXlS/9qj6e5rnNxUzcACFiVlITHtEa5bXi+mAmQbMwAFskHWGHHN2SOzmncExXNaZlZm0cnZt3XhCr/fhd87rwePLtT3DpCgEOIB18uxojDk7JKZfe3gMK6trF7+XNJ6k11t4y+W1PcOkKHShAFg3291kllhDnp/r6Njintyv1W/5G/X7Q88pYzxNQBcKgETDNemyFwBjNu/EvNa0LUhmQYADLZdUk05S1AJgzOadmNeatgXJLKiBAy0XM5MtegFw3BGtSYuTZY6nqQhwoOVC14bNmOmCe6ULgMOdJ5fObtLK2bVgF0rbEeBAy4Va8aredp7UedKZndF9+3YS1gFja+BmdrWZHTWzU2Z20szu6D1+j5ktm9mJ3n/vK3+4AIpW9U3qIW08zzuvmBn4OUl3ufvTZvZGScfN7HDve/e5++fLGx6AKtTh2rA2nued19gZuLufdvene5//RNIpSfx7BkCh2nied16p2gjNbLukXZKe6j30CTN71sy+bGabAz+z38yWzGzpzJkz+UYLoDb653Bfs/h4IedtF3m5Q1tE78Q0szdI+hdJf+juD5nZlZJ+qO7pjX8gaau7f2zU72AnJlC9Mm5qTzqHu4iFzzLGOg1COzGjAtzMZiU9JukJd/9Cwve3S3rM3a8b9XsIcKBaoy48mM8RkDHb4ZPGQjhnEwrwsYuYZmaSviTp1GB4m9lWdz/d+/KDkr5T1GABFCN04YGU7pCqYWkXHPMckIWwmBr4bkkflbRnqGXwj8zs22b2rKQbJN1Z5kABpDeug2N17bwOPHAidQ077YIjLYLlGDsDd/dvSBvuCpWkrxU/HABFCu2yHJZ2Rhxz7vcgWgTLwWFWQE0V0eWR1NkREjMj7o/pzgdOjL1xfhAtguVgKz1QQ0XVjAcvPFheWd1wU/uwUTPiPFvd087YEYcAB2poXM04TTfH4C7LfidIqKwyakYcexVaaAxpx43xCHCgQGXftt6fiWedmffDPNTHPWpGnLeOXYft+tOGGjhQkKpuWy+imyPLAVbUseuHGThQkDwlhmGhmnHogoMs3RxpZ8TUseuHAAcKUmSrXKhmHKpfXzXXKX2nI3Xs+iHAgYKEeq6zlhhCi4/DnSSd2RndcO2WSnY6UseuF2rgQEHSnKaXpsd7sLYudcO7v7OuX7s++twZdjq2EDNwIKfQPY6hEkPaHu/QeSbzY8oqEjsdpx0BDuSQZXNL2sXO2JbCJEV0iHCKYH1RQgFyyHJIU9rFzjQthYOK6BApsjUSxSPAgRyydJ6M6qdOqo2HauvnR5zlX9TFxJwiWG8EOJBDls0toUDud5IMz3YlJW66mQ+8xmBtPO91Z5wiWG/UwIEcsmxuGdXjHZrtHlvckzibTnrtIlsKi26NRLEIcCCHrJtbkvqp73zgROJzQ7PdLG8ERe0IZfdlPRDgQE5FbW7JMtst4o1gFHZf1hsBDtREUbPdMneEol5YxARqIssJgUnS7AhFszEDB2qkiNkuZY/2IMCBEZq6C5GyRzsQ4EBA3nspmxr+aA5q4EBAnl2IbEFHFQhwICDPLkS2oKMKBDgQkOcOSLagowoEOBCQpx2PC4BRBRYx0SixC4NFLCDmacdjCzqqQICjMWK7QvJ2jwzK2o5HLzaqQICjMWIPaSryMKc86MVG2QhwNEbswmCWBUR6ttFEBDgaI/aQprSHORVZcgnhDQJloAsFjRHbFZL0PFM3mJNup4nt2U667iwGm3pQFmbgaIzYhcHB5y2vrMok9W+PTJpdx5Rc8szSx9XkmZ0jKwIcjRK7MNh/3u6DRzaUUwZn1/c+8bxCVwMPllzyLIyOeoOoonyD6TW2hGJmV5vZUTM7ZWYnzeyO3uOXm9lhM3uh93Fz+cMF0gmFZz8ok2rl0sbSTJ6dlaM29bDlHnnE1MDPSbrL3d8m6d2SPm5mb5e0KOlJd3+rpCd7XwMTk1SjDoXnjNmG4OxLukghz87KUbV7ttwjj7EB7u6n3f3p3uc/kXRK0ryk2yTd33va/ZL2ljRGYKzQQuEN125JDM/znlw4MSnxBvi0C6ODRt20w5Z75JGqBm5m2yXtkvSUpCvd/bTUDXkze1PgZ/ZL2i9J27ZtyzVYICRUijj63Bl97kPXJ97cnqbVMO3CaNLPJ32PLffIIzrAzewNkh6UdMDdf2xmUT/n7ockHZKkhYWF0HoRkMuoUkQoPNMGZ8zCaJXnrQBRAW5ms+qG91fc/aHewy+b2dbe7HurpFfKGiQwTtrNO6OCc1xbX9F1a7bcI6uxAW7dqfaXJJ1y9y8MfOtRSbdLOtj7+EgpIwQiZClFJAVnTFtf2jcLoCwxXSi7JX1U0h4zO9H7733qBvdNZvaCpJt6XwOFi9kBOWqhMI2Ytr4854QDRRo7A3f3b6i74J7kxmKHA6yXZqNLEaWImPIIdWvUBTsxUWtVHw0bWx6hbo064DAr1FK/bBLaKZl1wXBcOYbyCJqEGThqZ7hskiTLgmFMOYbyCJqEAEftJJVNBmWdEceWYyiPoCkIcNTOqPLIfI4ZMeeOYNoQ4Kid0ELi/FxHxxb3XPw67Tna9G9j2rCIiYlKWlSMWUjMcssNC5SYNgQ4JiYUwpLGbsrJco52UZt9gLqghIKJGRXCSUe6Dspaz2aBEtOEGTgmpqxbboC2IMARlPUW9lhl3XIDtAUBjkRZFgnTyhPC1LMBauAIqOIMkry7Hqlno+0I8JbKemnB8sqqrll8PNcW87T92wCSEeAtlOfSAkkbWv4Gw3dcOKc5HhbAaNTAWyjrpQXDhn8mpm6epX8bQDICvIViLy0YXCSM+V0x4cx5JEBxCPAWim3f27trXscW9+h7B9+v+YifiQln+reB4hDgLZSlfS/mZ2LCmf5toDgEeAtl6aGO+ZmkcDZ1a+H9jUD0bwPFMXev7MUWFhZ8aWmpstdD9fpdKMsrqzJ1O1b6OrMzhDWQgZkdd/eF4ceZgaNQ/br5/FxHw1OD1bXzOvDAiVK25QNtRB84SjGqqyRrDzmA9Qhw5JYUvKM2Akkbt+WzwQdIjxIKcglt3rnh2i1jNwKl7SEHsB4BjlxCwXv0uTMXu01C0vaQA1iPAEcuo4K3v6D5xX07C+khB7AeAY5M+pc9hJpQB4M3aw85G3yA0VjExLpFyMs6szKTVs6uBTtBhhcchyUF77izu/OeDQ60EQHecsNhvLK6dvF7oU6QpLp333yO4OWCBiAdSigtNyqMpeROkFDd26Sxt8kDKA4B3nIxXR7Dz2HBEagHArzlYkJ3+DksOAL1QIA3WL8T5JrFxzOfLzLu5p3QgiQnCgKTN3YR08y+LOlWSa+4+3W9x+6R9DuSzvSe9ml3/1pZg8RGsVvPx50vMtz9EdOF0v85AhuYrLHHyZrZL0v6qaS/Hgrwn7r759O8GMfJbpT1AKfdB48knjUyP9fRscU9F3/3cLsfR7oCzZP5OFl3/7qkH5UyqpaLuQQ4JGbrOeeLANMtTx/4J8zstyQtSbrL3V9NepKZ7Ze0X5K2bduW4+Wmz7iAHTUzD532x/kiQHtkXcT8M0m/IGmnpNOS/jj0RHc/5O4L7r6wZcuWjC83nUJB2p+Jj5qZF3VHJYDmyhTg7v6yu5939wuS/kLSu4odVjuEgnTGbGzpg/NFAGQqoZjZVnc/3fvyg5K+U9yQ2uPum3ckLjKGdkYOz9g5XwRot5g2wr+V9B5JV5jZi5I+K+k9ZrZT3Ttrvy/pd8sb4vQKBWz/UuBhV811Unet0O4HTC9upa+hUPvfh985rwePLye2BUrjZ9rcOQk0U6iNkACvqaQjXl89u5b43LnOrF47d2FdsJu6/zzqnw4oiZ5woKEI8IYad/Z2jM7sjC6d3ZT4BjC48QdAPYUCnPPAa27cca8xVtfORy+MAmgODrOquXEB25md0ebXz2b+/fSEA81FgNfcqIDt935/9tffMfJEQalbJ6cnHJgulFBqLtQrnrT42G8/7C9gDj7/ng+84+Jz6EIBpgOLmA2Qtv2PdkFgutCFUiACEkCV6EIpSOxFCgBQNgI8pTxHwAJAkQjwlMYdAZtnZk5pBkAatBEOGXdRcJ4jYMe9btbbeQC0EwE+ICZEQ2dsnw8sBr+0shp1ezzXnwFIiwAfEBOioYsU5gMz88s6s1Eza64/A5AWNfABsSEaOmM7acONmYJvCmnvuASAQczAB2S5Q7JfHrnzgRN63SWbtPn1s+tm5iuBI2CH3xS4/gxAWszAB4S2rd99847EDhFp/ax7ZXVNndkZ3bdv57rbdmJm1lx/BiAtdmIOiQlqKf6M7dDtOlykACBWa7fSF9FbvfvgkcRZdIhJ+t7B9yeOoX+7zsrZNWbZAKKEAnyqa+BF9Van7QRJKo8cW9yj+/bt1GvnLujVs2v0egPIbapr4KG2wAMPnNA9j568OBMeNysOdYgk3UU5auFxVJsis3AAaU11gI+aOa+sriV+PrgFXir2jG16vQEUaeoCfLDevMksuENylNW187rn0ZPrZteujTe994M6dvZMrzeAIk1VgA93fGQJ777BWXlfP7xjbnFPWjwd1aYIAGlN1SJm6Ab3GbPCXiOm3BFaPJWUuA2f+jeALKZqBh4K1wvu+uK+nRtmvyGjerxjyh2jFiuPLe4hsAEUYioCvF+uCBVMrprrbNjpONh5ktSFIiVv3okpd7BYCaAKjQ/wpJ2OgwZDN3QI1ShZNgGxWAmgCo0P8FDdW9rYLZJWlsCXRp+pAgBFaXyAh8oSJkV1i5SBg6kAVKHxAV7XckXW2TsAxGp8GyHnaANoq8bPwClXAGirxge4RLkCQDuNLaGY2ZfN7BUz+87AY5eb2WEze6H3cXO5wwQADIupgf+VpFuGHluU9KS7v1XSk72vAQAVGhvg7v51ST8aevg2Sff3Pr9f0t5ihwUAGCdrDfxKdz8tSe5+2szeFHqime2XtF+Stm3blvHluoq4Hg0ApkXpi5jufkjSIal7J2ban++H9vClCoMn/BHiANooax/4y2a2VZJ6H18pbkg/M3gsq6QNh1X1T/gDgDbKGuCPSrq99/ntkh4pZjjrjTrnpG95ZVW7Dx7hYmAArRPTRvi3kv5V0g4ze9HMflvSQUk3mdkLkm7qfV242ONXud0dQBuNrYG7+0cC37qx4LFsEDrnJAm3uwNom1qfhZJ0zsmoy9G4MAFAm9Q6wPfumt9wh+R9+3ZqPnDS4KRPIASAKtX+LJTQOSdcmACg7Wof4Ek4gRAAGhrgEicQAkCta+AAgDACHAAaigAHgIYiwAGgoQhwAGgoc099wmv2FzM7I+kHlb1gca6Q9MNJD6JCbft7Jf7mtmjq3/wWd98y/GClAd5UZrbk7guTHkdV2vb3SvzNbTFtfzMlFABoKAIcABqKAI9zaNIDqFjb/l6Jv7ktpupvpgYOAA3FDBwAGooAB4CGIsBTMLNPmpmb2RWTHkvZzOxeM3vOzJ41s6+a2dykx1QWM7vFzJ43s++a2eKkx1M2M7vazI6a2SkzO2lmd0x6TFUwsxkze8bMHpv0WIpCgEcys6vVvcD5vyY9looclnSdu/+ipP+Q9KkJj6cUZjYj6U8l/Zqkt0v6iJm9fbKjKt05SXe5+9skvVvSx1vwN0vSHZJOTXoQRSLA490n6fcktWLV193/0d3P9b78N0lvnuR4SvQuSd919/909/+T9HeSbpvwmErl7qfd/ene5z9RN9Sm+nB9M3uzpPdL+stJj6VIBHgEM/uApGV3/9akxzIhH5P0D5MeREnmJf33wNcvasrDbJCZbZe0S9JTEx5K2b6o7gTswoTHUajG3shTNDP7J0k/n/Ctz0j6tKRfrXZE5Rv1N7v7I73nfEbdf3J/pcqxVcgSHmvFv7LM7A2SHpR0wN1/POnxlMXMbpX0irsfN7P3THg4hSLAe9z9vUmPm9n1kq6R9C0zk7qlhKfN7F3u/j8VDrFwob+5z8xul3SrpBt9ejcMvCjp6oGv3yzppQmNpTJmNqtueH/F3R+a9HhKtlvSB8zsfZIulfRzZvY37v6bEx5XbmzkScnMvi9pwd2beKJZNDO7RdIXJP2Ku5+Z9HjKYmaXqLtIe6OkZUnflPQb7n5yogMrkXVnIvdL+pG7H5jwcCrVm4F/0t1vnfBQCkENHCF/IumNkg6b2Qkz+/NJD6gMvYXaT0h6Qt3FvL+f5vDu2S3po5L29P7fnujNTtEwzMABoKGYgQNAQxHgANBQBDgANBQBDgANRYADQEMR4ADQUAQ4ADTU/wMeoX4vTceKnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make sythetic data (from Eric's example)\n",
    "x_data = numpy.linspace(-5, 5, 100)\n",
    "w_true = 2\n",
    "b_true = 20\n",
    "\n",
    "y_data = w_true*x_data + b_true + numpy.random.normal(size=len(x_data))\n",
    "\n",
    "pyplot.scatter(x_data,y_data);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This situation arises often. In **Module 1** of _Engineering Computations_, we used a real data set of Earth temperature over time and we fit an ominously sloped line. \n",
    "We derived analytical formulas for the model coefficients and wrote our own custom functions, and we also learned that NumPy has a built-in function that will do it for us: `numpy.polyfit(x, y, 1)` will return the two parameters $w, b$ for the line\n",
    "\n",
    "$$y = w x + b $$\n",
    "\n",
    "Here, we will instead use gradient descent to get the parameters of the linear model. \n",
    "The first step is to define a function that represents the _deviation_ of the data from the model. \n",
    "For linear regression, we use the sum (or the mean) of the square _errors_: the differences between each data point and the predicted value from the linear model (also called _residuals_).\n",
    "\n",
    "<img src=\"../images/residuals.png\" style=\"width: 400px;\"/> \n",
    "\n",
    "#### Each data point deviates from the linear regression: we aim to minimize the sum of squares of the residuals.\n",
    "\n",
    "\n",
    "Let's review our ingredients:\n",
    "\n",
    "1. observational data, in the form of two arrays: $x, y$\n",
    "2. our linear model: $y = wx + b$\n",
    "3. a function that measures the discrepancy between the data and the fitting line: $\\frac{1}{N}\\sum (y_i - f(x_i))^2$\n",
    "\n",
    "The last item is called a \"loss function\" (also sometimes \"cost function\"). Our method will be to step down the slope of the loss function, to find its minimum.\n",
    "\n",
    "As a first approach, let's again use SymPy, which can compute derivatives for us. Below, we define the loss function for a single data point, and make Python functions with its derivatives with respect to the model parameters. \n",
    "We will call these functions in a sequence of steps that start at an initial guess for the parameters (we choose zero), and step in the negative gradient multiplied by a step size (we choose $0.01$). \n",
    "After $1000$ steps, you see that the values of $w$ and $b$ are quite close to the true values from our synthetic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHgAAAAaCAYAAAB8WJiDAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE+UlEQVRoBe2a7VHcMBCGjxsKIKSCQAcQKgjpIEAFQAcw/OMfAx0AFRDoAKiAjw5IByHXAXkfI+1Isuz4jvMRJ94ZnaTVStoP7WplmHt5eRn00G0NHB4eLkiCAyfFkqu3hR/Nd1u0nnungWMZc9drQ+1TtR9Uloce2ded1sCOjLoeSHCs9pJwK72BA610uIn33uf4n+vv4Jxauo2T5+LB31SXQ7SQ/pLutpT/EPfj2ES0KxhXZRUVRCFag3vCQdDD36UB7lNsUwvuIOC9q2qPILYQLQRWX1O9zwCgNun3pcpnlWf1l1X38A4akO4x8Ej1WW574Ym8+6qLbNr1Xw2sDoa8VV24dbqA8KTc96otFU9p+n77GnB2+KK68E6/o/oYl6eROafa2Grfh2jcGoIqIGxfVw02xYuRdZWdpvQ9XUkD2AhbpYAD8kyi9oWnk33o2FQn650Yxa124+q3VEQKSg8TaEC2OFP5pUIoNi9W+0PVckMNcvf+qCIQ/ivj4YI1tP1Q+xrAVptNt+FTJQas8048+FEG9lncmvoX6l+pngloL7ye0P5R5Ul9SzTU5vrYVe2TCw5slEmqP9C4X4MuMmyrcHdtqQB3onkPmUhcH7R3KBO8ngu3AWMJYCtsZvTJeNTlDiZDfoqwruOUggJRBOHhRDWKuVQb3KzgwO19pw3TO+hAuJAXQhd95AqB77UnwTrnGiQnIDFBoem64dw22l4mcpt0bzyUg5oDbBXKm6Mx3FAthHs2TNzw969lblIICqRUMRCv8Mae9uOAYViAk5vyCo+WAIqeE/6oYteOcHh/qETPv/eCxWRc3faggUzIWRVVkb+xgQnRCIfAOWAjwnM6viA84TILoifb84cjpGGvgcZzCR375EIS9z8GAzjZR0VLP8JjfHhJlcGTzgys8bTPc9DkEm1uX5G8gsbZ41aFuilsaJ7nO51TKZMjRHcmZzIZuRrzgYHrgI2ie0lMe8N5ryrNF03OgAPh8Xq+yhDqG4FoRxC6uQjmvQ40vPD4TxUZXTmZ8eigsFAdaD48ZL8R1M2rGnPrDVTDfySTcFWH1i9X55CexmpCNC7PJjkgFKSG9Kc99Zrc/GnictGkFMqcglKDGx+BUu3gCrdAMaLZNdClRRK3bdWh9VzBJzZrBPOiwuXrYrqFOqcETj9Z66jRDtMjgkfjxS0LzgzlcFvizb7oOJ4vNUaSxaEkuuD14VokPDbHrTOLKidT6dAmjODBIe/JcNwdqstpX4vR1kN5YTaKor5LGWGYNOKWG/BpXiYeOOkIi5IKcDhLuBwaOsqzxkun381Jo5Sb2noVGcrxAq+pDCEjRf4QIura8xq8UMFwOdgWkvcYbzWUw/PoPYw70L58vTmliA+ix0+VTyrw53HRG1ljAF4LzyhuIFqiz7Wbw2c9DJ9GAUhnAUQN+CfDRx7/x5y66w85sjmO8CUo/pqkDUhK6rK+0sRJENpn7CRrkn26OscZel11NqETnmjF4fQH4Y+iEqIBTlDjU1HMmOwHz6P894AxVYggBahNhOS9fvSKyf7i8eF7PksUIgsDa3FCGM8Xu89Comm1tf6N22taS3Z5HaIZ16OHczX4Wpi9LpxtsNFYVyR3sAdSdu5isrge2tcA3rgog+2pJuSSX9TdveQZY0dZ+48OTR5oAzyYf9Zq/CGCeT20qwHZg0NwpTrKupvs+ht4KM5n5o25AwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle \\left(b + w x - y\\right)^{2}$"
      ],
      "text/plain": [
       "             2\n",
       "(b + w⋅x - y) "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w, b, x, y = sympy.symbols('w b x y')\n",
    "\n",
    "loss = (w*x + b - y)**2\n",
    "loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "grad_b = sympy.lambdify([w,b,x,y], loss.diff(b), 'numpy')\n",
    "grad_w = sympy.lambdify([w,b,x,y], loss.diff(w), 'numpy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be sure to read the documentation for [`sympy.lambdify()`](https://docs.sympy.org/latest/modules/utilities/lambdify.html), which explains the argument list.\n",
    "Now, we step down the slope. \n",
    "Note that we first compute the derivatives with respect to both parameters _at all the data points_ (thanks to NumPy array operations), and we take the average. \n",
    "Then we step both parameters (starting from an initial guess of zero)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.03136593978268\n",
      "20.030956250406764\n"
     ]
    }
   ],
   "source": [
    "w = 0\n",
    "b = 0\n",
    "\n",
    "for i in range(1000):\n",
    "    descent_b = numpy.sum(grad_b(w,b,x_data,y_data))/len(x_data)\n",
    "    descent_w = numpy.sum(grad_w(w,b,x_data,y_data))/len(x_data)\n",
    "    w = w - descent_w*0.01 # with 0.01 the step size\n",
    "    b = b - descent_b*0.01 \n",
    "\n",
    "print(w)\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAobUlEQVR4nO3deZzNZf/H8dfHGJpKlpKYiCQUoqR+TcrS5m4xuEu7FtTdRiRLi9J9Z0riV9x1U35p0V2JCVkSUrQSFdnSPoSWEZms1++PMzPNcr5nm3Nm5sx5Px+PHsyZc77fax49Hu+5fL6f67rMOYeIiMSfSmU9ABERiYwCXEQkTinARUTilAJcRCROKcBFROJU5dK82RFHHOEaNmxYmrcUEYl7y5cv/9k5V7vo66Ua4A0bNmTZsmWleUsRkbhnZt/5e10lFBGROKUAFxGJUwpwEZE4pQAXEYlTCnARkThVql0oIiKJJnNFFqPmrWNTdg71aqQw6PympLdJjcq1FeAiIjGSuSKLodO+IGfvfgCysnMYOu0LgKiEuEooIiIxMmreuvzwzpOzdz+j5q2LyvUV4CIiMbIpOyes18OlABcRiZF6NVLCej1cCnARkRgZdH5TUpKTCr2WkpzEoPObRuX6eogpIhIjeQ8q1YUiIhKH0tukRi2wi1IJRUQkTinARUTilAJcRCROKcBFROKUAlxEJE4pwEVE4pQCXEQkTqkPXESE2G77GisKcBFJeLHe9jVWVEIRkYQX621fY0UBLiIJL9bbvsaKAlxEEl6st32NlaABbmYHmdnHZvaZma02swdzX69lZvPNbEPunzVjP1wRkeiL6bavW7ZA//6wfXvJr1VEKDPw3UAn59xJQGvgAjM7HRgCLHDONQEW5H4tIhJ30tukMrJ7S1JrpGBAjZRkDkquxJ2vrCQtYyGZK7LCv+ju3fDoo9CkCYwfD+++G/VxBw1w57Mz98vk3P8c0BWYnPv6ZCA96qMTESkl6W1SWTqkE2N6tmb3vgP8tmsvjr86UkIOcecgMxNOPBEGD4azz4bVq+Hii6M+5pBq4GaWZGYrga3AfOfcR0Ad59xm33jdZuBIj8/2NbNlZrZs27ZtURq2iEhslKgj5fPP4ZxzoFs3qFoV5s2DmTPh+ONjMtaQ+sCdc/uB1mZWA5huZi1CvYFzbgIwAaBt27YukkGKiJSWUDtSCi78OaHybsatnU6j6VOgRg0YNw5uugkqx3apTVhXd85lm9k7wAXAFjOr65zbbGZ18c3ORUTKtYLBWz0lGTPI3rU3f/VlvRopZPkJ8YIdKXkLf/b9+Sc3LJ9Fv/f/y8F7ctjYsxeNxz8GtWqVys8SShdK7dyZN2aWApwDrAVmAL1y39YLeCNGYxQRiYq84M3KzsEB2Tl7i9W6OzarHbQjZdTctfzPmg+YO+k27lv0LJ/Wa8b5N4yn8zE9SJuwMrKHnhEIZQZeF5hsZkn4Av9V59wsM/sAeNXMbgS+By6N4ThFRErMX327oJy9+1m0dhsju7f03hdl9WoyJgyi/bcr+KrW0Vz39+G80/jU/GuU5jL8oAHunPscaOPn9V+AzrEYlIhILISysnJTdo7/g4h/+YWvbxlIg6kv0Co5hQc79+GFNheyL6l4jOY99CzzABcRqSi86ttF31PI3r3w73+z5/7hNNixg5fadGHMmVeRnXJYwOuUxjJ8LaUXkYThb8VlQcVWX86eDS1bQv/+rDzyOLpc/yTDz/1H0PCG0lmGrxm4iCSMvJJGoC6U9DapsGYNDBgAc+f6VlLOnEnP98CZ+b1uSnJSodp61JbhB6EAF5GE4re+nefXX6FfP9/S90MOgcceg9tvhypVqLdqod/yS2pu8JfFYRAKcBGRffvg6adh+HDIzoY+feChh6B27fy3DDq/aaFDH+CvmXbAXwoxpAAXkcQ2fz7ceadvv5KOHWHsWGjVqtjbipZfysOxawpwEUlM69fDwIEwaxY0bgzTp0PXruBR54Yg5ZdcpXm2pgJcRCq8gqF6fNX9jN84i+Ne+T9ISfFt+XrHHb7Np6Jwn9I8W1MBLiJxJdwZbl6o7tm9hys/m8eA916kZs4Ovu12OQ2fGgN16kRtbIF2MlSAi0hCi2SGO2reOtpsWM59C5+h+bZv+bB+C0Z07sv2pieyNITwDucXRmmfrakAF5G4EfYM96uvGD7pHs7b8CE/VK/DzelDmXv8GWCGhRCq4f7CCGUnw2jSSkwRiRshz3C3b4e774YTTuDM71by6FnXck7vp5jbNC3/IWUooRru4Q4xPVvTD83ARSRuBJ3h7t8PkybBvffC1q1w3XUsvuJ2/u+9reyOYKVkuCWR0m41VICLSNwItJiGxYt9p7+vXAlpab59TE45hS7A7tqRtfZFUhIpzUU9CnARiRv+ZrjDW6Rw3r/6weuvQ4MG8N//wmWXFernjjRUA/7CKAcU4CISM6F0cITbFpgfxjt2wMMPQ4/HfWdPPvSQb2FOSvQeGJbH1ZcFKcBFJCZC6eCIaOHLgQMweTIMGwY//QTXXusL8tTUmKyCLKt9TkKhABeREvMXnKG0/IXdFrhkiW+3wE8/hdNPhzfegHbt8sdQmqsgywO1EYpIiRQ9KDgvOL1OvinYwRFyl8d338Hll0P79r7ukpdegvffzw9vCL/lryJQgItIiXgFZ5LHplAFOzi8ujnyX//jD7jvPmjWDGbM8G33unYtXHllsU2nSnsVZHmgABeREvEKyP3OBV3U4rnw5dwm8PzzcPzx8M9/QvfusG4dPPCA76AFP4L+MqiAFOAiUiKBArJq5UrUPDgZw3dyzcjuLQvVo9PbpDKye0tSa6Tkv+fpxntIv+Xv0KsXpKbC0qW+kkn9+n7vkbkii7QM32k5Ref8hq+kk5axkMwVWSX+Wcsbc86V2s3atm3rli1bVmr3E5HYK/rwsKiU5KRiwe3XDz/AkCEwZQrUrQsZGXD11VDJe57p794GuAJ/hj2OcsjMljvn2hZ9XTNwESmRgrNof4I+SNy1y1caadoUpk3zLYNfv97XHhggvMF//d0BSWYUnZpWxAeaaiMUkRLL65VuNOTNYsEJxevkmSuyGDV3LW0/mMuwdydTZ/s23+rJRx6Bhg1Dvm+g+ns4749XCnARiZpQ9g7JXJHFlCen8uS8pzh50zq+qNOYgV3v5u/9ryC9YXjlDa/7JZn5DfGK9kBTAS4ixUS6ojHo3iFZWVS98XpeXTGfrYfUZFCXfkxt2RlnlfimyOKdgmOonpKMGWTv2ltoPF7363FKKq8vzyq3e5hEiwJcRAopyYpGz71DmtXytQOOHEmn3Xv49+l/Z/zpl/FH1YPzP1uwvFF0DNk5e/O/5288/n7ZtD2mVrndwyRa1IUiIoVmu5U8yg+pNVJYOqRTeBd2Dl591Xe4wvffQ48eXHpsOp9Uqhnw+nltgYFENJ44pS4UEfGr6FL4qD0AXL4czjrLtwS+Zk1YtAimTuWqKzoGXeATyr0q2gPJSCjARRKcv1Y8f0J+ALh5M1x/PZx6qq8dcOJEX5h36AD4X7xTtD87lHtVtAeSkVANXCTBhTKTDekB4J9/wpgxvq1dd++Gu+6Ce+6B6tWLvTXYFq3+Hk6GPZ4EoAAXSXCBWvEOOBf8AaBzvtNwBg2Cb7+Frl3hscfguOPCHkvRzpODkiuRvWuvZxdKolOAiyQ4r1a8kJadr1zpO4dy8WJo2RLefhs6d45oHP46T1KSkxjTs7XC2kPQGriZ1TezRWa2xsxWm1m/3NcfMLMsM1uZ+9/fYj9cEYm2UGrSxWzZAn36wMknw+rV8NRTvkMWIgxvSMz9vEsqlBn4PmCgc+5TM6sGLDez+bnfG+Oceyx2wxOR0hDysWG7d8MTT/jOn8zJ8c2+778fatQo8RgScT/vkgoa4M65zcDm3L/vMLM1gP49I5JInPMdX3bXXbBxI1x4IYwe7duAKkpCWYYvhYXVRmhmDYE2wEe5L91mZp+b2SQzK96Z7/tMXzNbZmbLtm3bVrLRikjp+/xzOOcc6NYNqlaFefNg1iwydx1KWsZCGg15Myr7bXse7qBuE08hr8Q0s0OBxcC/nHPTzKwO8DO+3RsfAuo6524IdA2txBQpfRGf1L5tm+84s4kTfSWSESPgppugcmW/+3BHY7/tWJwqXxF4rcQMKcDNLBmYBcxzzj3u5/sNgVnOuRaBrqMAFyldgQ48SPUKyD17YNw4X2Dv3Am33OLbr7tWrfy3eC11D7S8XeEcOa8AD1oDNzMDngXWFAxvM6ubWx8H6AasitZgRSQ6vA48AD+bQjkHs2bBwIGwYQNccAE8/jg0b17suuE+cCzJBlniLZQaeBpwDdCpSMvgo2b2hZl9DnQE7ozlQEUkfME6OHL27qf/Kyu5uv+zbD2jA1xyCSQlwezZMGeO3/CG8A8QVotgbITShbIEip0VCjA7+sMRkWjy6uzIU3PXdu5cMoWrVs5hZ9WD+fyuB2j18DBITg543aD7fhehFsHY0EpMkXIqGjVjrz1FKu/fxzUr3qT/kikcuieHl1p3YcyZV3Lw4UexNEB4ey11DzY+tQjGhgJcpByKVs244IEHWdk5GHD2xk+4b+GzNP71R947pjUPde7N+toNAcgOMCMuyVL3cGfsEhoFuEg5FKxmHM7MPH+V5Zo1/NTnVo5auoiva9bjxh73saBxO7C/KqSBZsSBxhTxST16gFkiCnCRKIpWq5xXbThvJh7WzPzXX+HBB2H8eI469FAYPZovzuzG+zPXQRgz4pLWsUNeri8h04EOIlFS9GSbvHCNZIWi10w4ySz0bo59+3z93E2a+P7s3dvXHjhgAF3bNQp7A6twO08k9jQDF4mSkpQYivKqGXsdcFBsFvzWW3DnnfDll9Cxo++ghZNOKvSWcGfEqmOXP5qBi0RJNFvlvLZ4TQ0wC85ckcXlAyez4Lh2cP757Pz9D5g+HRYsKBbekYho21mJKc3ARaIk2q1yBWfIebX1vE6SghtgpCQn0aV+VbL/cTsvfDKDPytXYWSH63j59G6MOOYU0s3fMo7IqI5dvijARaIknBJDOA87i7bvOf7az6RBtSo8sWs5x9yUQfU/fueVVucy+qxr+PmQmuCIqHwj8UMBLlJC4S5uCbfH22s/k4t/Wcs/X3+G6l+t5aP6LRhxaR9W12lc6H1a6VixKcBFSiCSxS3hPuwsGsLH/LaJYYsmcf6GD/mxeh2GdB3CnKZphfq580SjQ0S7CJZfCnCREoik8yTch515tfVqu//g1vdf4YZlM9ibVJlRZ/fimbZd2V25it/PRaNDRLsIlm8KcJESiKTzJNDDTn+z3UHnHMenw0dz+zvPU3tXNq+1OIcnO1/P9wdV97yH517fYYpma6REn9oIRUogksUtXkeHdWxWu9hCoGljXqLD1X9jxOwn2FT7aC65dgxjrxrKgF5ne7YU5oX3qHnrSnzcmXYRLN80AxcpgUgWt3jtC1Jwtnt09k8MWzSJv61/n59qHAn//S8nXXYZM4rUuf3dO+8XQTTKHtpFsHxTgIuUQKSbNPnrp77zlZUcsnsXt3z4Gr0/yWR/pUo81v5qnjm1G2t7dg/53qWxIlSrL8sHBbhICUVlccuBA/T+ajF95j7DkX/8xusnduTRs3uxpdoRnqUSr3vf+cpKv++NdEUoaBfB8koBLlLWliyBfv2459NPWZnajL7d72VlPd8MN5LZbixXhEr5ooeYImXlu++gZ09o3x62boWXXuLbGfPZdkLrEu014vWQVGWPikczcJHStnMnPPIIPPaYb/HN8OEwaBAccgjpQPrJR5fo8ip7JA4FuEgAUV2FeOAAvPACDB0KmzfDlVdCRgbUrx/dQaOyR6JQgIt4KOkqxILhf972r8l471lqrv4MTj0Vpk6FM86I6fil4lMNXMRDsHMpA8kL/wPff8+YmaP4z9N3sOf7H1k+Yix8+KHCW6JCM3ARDyVZhfjEzM+4adGL3PTRNAzHE//Tk6dP/zs1k2uxtJLmTRIdCnARDxG14zkHL7/Mi4/1o96On5nVrD0jO1xPVvUjAcjREnSJIk0FRDyE3Y738ce+0shVV/F7tZpcemUGt3UdnB/eoCXoEl2agUtcCbUrJBrdIyG342Vl+TpLXngBjjoKJk1ibatzWJW5GrQEXWJIAS5xI9SukGjuYR2wHS8nx9fLnZEB+/f7QnzoUKhWjXSASpXUiy0xpQCXuBHqJk0x38PaOXj1Vbj7bvj+e+jRAx59FI49ttDb1IstsaYAl7gRaldIJN0jIZdcli+H/v19+5ecdBJMngwdOoT6I4hElQJc4kaoXSHhdo+EVHLZvBnuuQeeew6OOAImTIAbboCkJL/X9HcPlVMk2tSFInEj1K4Qf+8zfMHs73SagAt2/vwTRo6E44/nwAsv8lL7S2l11TjSfmlM5uc/hTTuvF8QBU/aGTrti4hPyRHJoxm4xI1Qu0IKvi8rOwcDXO73/M2u/ZZWnKPVR2/DCb3hm2/Y1PF8rm9+Geuq1QHg9zAejAaryWt2LpFSgEtcCfXBYN770jIWFiunFFwOP2reuvxwz3Pilo3cv2Aip/2wClq2hAULuPRjPK8TbDyBavI69V1KImgJxczqm9kiM1tjZqvNrF/u67XMbL6Zbcj9s2bshysSHq/wzAvKgqF8xB+/MXLOE8x8rj9NfvmelcMyYMUK6NSpRMvqAx18XJL9VkRCqYHvAwY655oDpwO3mtkJwBBggXOuCbAg92uRMpO5Iou0jIWFTmL3Cs8ks/zgrLJvLzd9NJVFE/ry91ULeDWtOx/M+YDW/xqc/5AyktPn8wSq3evUdymJoAHunNvsnPs09+87gDVAKtAVmJz7tsngW7sgUha8HhR2bFbbb3judw6c47z1H/DWs7cw9J3n+Kh+Cy64YTyXL5nKhWedUOgz4T4YLSi9TSoju7cktUZKsZN2SvKLQSSsGriZNQTaAB8BdZxzm8EX8mZ2pMdn+gJ9ARo0aFCiwYp48SpFLFq7jZHdWxZ7SPj65DncnPkkad99zvrDG3DNZSN4r9HJngcIh/tg1N/n/X1Pp75LSZhzRR/heLzR7FBgMfAv59w0M8t2ztUo8P3fnHMB6+Bt27Z1y5YtK8l4RfxqNOTNYg8jwTdL/ibjwr9e2LYN7rsPN3Ei26sewuNpVzKldRf2JVUmJTkppDMo/T0YBd/MeumQTmGPXV0oEoyZLXfOtS36ekgzcDNLBl4HXnLOTct9eYuZ1c2dfdcFtkZvuCLhCbp4Z88eGDcORoyAnTux225jabc+LPhwK/uzc0gtEJzBAjXadWstuZdIBQ1wMzPgWWCNc+7xAt+aAfQCMnL/fCMmIxQJgWcp4rzjYeZMGDgQNmyALl1g9Gho3pwLgQs7FL5OKG19Ee0TLhIDoXShpAHXAJ3MbGXuf3/DF9znmtkG4Nzcr0Wizl93SVH+HhSOa5VM+pAb4JJLfN0ks2f7/mve3PNeobT1hb1PuEiMBJ2BO+eW4Csl+tM5usMRKSychS75pYiff4bhw+He/0C1ajB2LNxyCyQnB71fKOWRkPcJF4kxrcSUci2srWH37oV//xseeAB27ICbb4YHH4TDDw/5fqGWR1S3lvJAm1lJuZRXNvEXpuBnpjx7tm/Ze//+0K4dfPaZ76FlkfAOVo5ReUTiiWbgUu4ULZv4kz8jXrMGBgyAuXOhSROYMQMuugiseNUvlHKMyiMSTxTgUu74K5sUlJKcxLDTj4R+/WD8eDj0UF9nyW23QZUqYV3XXzlG5RGJFwpwKXcC9VM3qFaFJ3//iJN6XAXZ2dC3r6+3u3btiK+rfUckXqkGLuWOVz91t62reHdKf0569D5o3ZqFU+aSdsylNBr9cdD9SAJdV/3bEq80A5cy5W/VY9FFOcf+8iP3vzOJDl99DI0bQ2YmmfVPYej0VeTs3QOEth+J9h2RikYzcCkzXjsIAozs3pKmVfdz74KJvDXpVtI2f+k7+X31aujalVFvrQ97H+1AuwKKxCPNwKXMeD1UfHz2l7xb6yvSx98Hv/4KN94I//wn1KmT/75I69l6QCkViQJcyoy/sP2f7z5j+NsT4Ofv4KyzfKso27Qp9j7tRyKiAJcAYr3NacEQPua3TdyzaBLnbfiQTTWPgqlToXt3v/3coHq2CKgGLh686tPBOj3CMej8ptTe/ydDFk1i/jO3cMZ3n/F4x+tYNnsJ9OjhGd6gerYIaAYuHsLagyQS+/eTvmw2F0waRpXffmFqi868eFFfbrj0DC4J8fqqZ0uiU4AnqEgPLcjKzqHRkDdLVFJZMvE16gwfSpPNG1nXsAW/PvE8l13Vhcsi/mlEEpMCPAGV5NACoFjLX8EQD/iL4euv2dTnNs5cOIcfD6vNrZcM5s1mZ5Ky1hi5IkuzaZEwqQaegCI9tKCoop/xqpvPWrIOhg6F5s2p+d4iHmt/NZ17P82bzduDWdD+bRHxTzPwBBTJoQVeR18X/EzRXwzmDnDR8vmc/r8vwI5f4Zpr6FitMz9VOyLkMYmINwV4Aork0AKvvbkLfqZgCJ/6wyruXzCRlls28mm9phwxfzacdhpJGQtB/dsiUaESSgKK5NCCUD5Tr0YKqdu3Mu6NR3htyhAO37WdOy6+i9tvGw+nnRbxvUXEP83AE1AkhxYE/czOnUzaOINjnnsKhzE27Qr+064Hf1Y5CLf9T9IyFhZ6vw5MECk5c86ruhl9bdu2dcuWLSu1+0kpOHAAXnzR95By0yZ+6JLOHa16soJqGBSqnackJ2mxjUgEzGy5c65t0ddVQpHIvf8+nH469OoFRx8N779P/dnTmZ5xOak1Uoo9+MzZu5/+r6wMae9uEQlOJRQJ3w8/wODB8PLLUK8ePP88XHUVVPprPhCoqySiHnIRKUYzcAndrl3wwAPQtClMnw733gvr1pHZohNpj75T6KT3YF0lofaQa6Yu4k0BLsE5B1Om+IL7wQfh4oth7Vp46CEyN2z3G7wdm9UOuhAoUA85BD+gQSTRKcAlsI8/hrQ0X4mkTh147z145RU45hjAO3gXrd2Wv1ugF68e8oK0wEfEmwJc/MvKgmuv9fVvf/01TJrkC/Mzzyz0tkDBm94mlaVDOjG2Z+uQesj90QIfEW8KcCksJ8d3fNnxx/tm2kOGwIYNcP31hR5SZq7IIi1joecS+4LBG8re3VrgIxI+daGIr/tj7lrafPgWwxY/R73tW5lz/BlM6voPrrqsI+nVqhV7f9HTcAryF7zB9u7WAh+R8CnAE1zmiixeGD+NMfOept2PX/LlkY0YcMXDfNigFQCr/LT7+at750ktQfDqgAaR8CjAE9nmzST3uZHXlr/FLwdXZ/AFt/Nay3M4UOmvUoa/U3i86t4GLB3SKdajFpFcCvBE9OefMGYMPPww5+76kwntujH+jJ7sqHqI37cXDWydCC9SPijAE4lz8PrrMGgQfPstdO3K1cf14OPKtQJ+rGgw60R4kfJBAR7Hwlp6vmIF3HknLF4MLVrA229D585cuSKLLyJ4IAl64ChS1oIGuJlNAi4CtjrnWuS+9gDQB9iW+7ZhzrnZsRqkFBfKuZYAc95eyd5h93DRJ3PYfshhfDdsJK0fvAsqVy703rwwrp6SjBlk79obMJj1wFGk7IUyA38OGAc8X+T1Mc65x6I+ogQT6QZOgZaep7dJhd27WXX3CNo/PYaq+/Ywqe0lPJF2BXsrVWfkF1sK3UNhLBKfgga4c+5dM2tYCmNJOKHOov3xXAH52y7IzIS77qLFxo3MP64dD3e8kW9q5V7PT1eJiMSnktTAbzOza4FlwEDn3G/+3mRmfYG+AA0aNCjB7SqeYBs4BZqZ++sEabb1G/61+Fl4dCWccALXXjaCdxudXOy+2l9EpGKIdCn9U0BjoDWwGRjt9Ubn3ATnXFvnXNvatWtHeLuKyStI82bigbZWLbj0vNau7fxr3jjefK4fLX/5FsaNg88+Y2ObNL/XV7ufSMUQUYA757Y45/Y75w4AE4F20R1WYvAK0iSzoFurprdJJePipgz4YhbvTOhLz8/f4tvLr6PK1xvh1luhcmXtLyJSwUUU4GZWt8CX3YBV0RlOYvEK2P0e55Tmz9idg5kz6Xr1+dwx+2kO69ieyqtW0XjKs1Drr57uUDaREpH4FUob4ctAB+AIM/sRGA50MLPW+M6s/Ra4KXZDrLi8+qlHzVvnudJxwWsLqDZsMO2+Ws63tRuw6cnnOeO2awLeQ4EtUjHpVPpyyN9uf3X37mTMl9M5dd5r7KySwtgzr+SFNheSfFBVRnZvCQRfWKMzJ0Xik9ep9ArwciovbLf+soM+q+Zy8zsvcvDuXbzUpgtjzryK7JTD8t9bIyWZ3fsOFAp8w/fPo7zdAQG/y99VUhEp/xTg8WjOHHbccgfVvv2Kdxu24aFOvdlQ+5iwL5OSnMRByZX4bdfeYt9LrZGiHQRFyjmvANdeKOXR2rUwYADMmcNvh6fSr8f9LGx8KphFdLmcvfs99zpRT7hI/FKAlye//uo79X38eDj0UBg9mnM2H8uepGTPjwSaXYdCPeEi8UtnYpYH+/b5Ft80aeL7s3dv3zmUAwZQ+/DDPD+W1xY4/OITi7UjFlUjJVk94SIVjGbgZe2tt3zbvH75JXTq5DtooVWr/G977b3t7+FjXvth3gPMgu9/4JIT89+jLhSRikEBXlbWr4eBA2HWLGjcGKZPh65di9W5Q917u2C/d6B2QQW2SMWhLpQIlKifOjsbRoyAJ5+ElBS47z644w6oWjWmYxaR+KUulCiJeAvYffvgmWd8gf3LL74690MPQZ06pTFsEamAFOBhimgL2AULoH9/WLUKzj4bxo6F1q1LfewiUrEowMMUbAvYgjPzp56Zx2lfvkrdd+ZBw4YwdSp07+7Zz62l7iISDgV4EcFC1N9BClB4C9hqu//g1vdf4YZlM9hXuTI8/LCv0+SggwLeN9LTeUQkMakPvIC8EA31IIU8eVvAVjqwn8tXzmXhhJvo+/F0Mk/sQIc+E8i84FrSxr5PoyFvkpaxsND18gQrzYiIFKUZeAFBDwrGu63v7ade4ZY3xnHC1m/4JPUErr/0AVYddRw1UpJDmll7nnGppe4i4kEBXkCoIVpoj+2vv4ZBd5A+bRqbDjuSWy8ZzJvNzgQzUpKTMCPoLwXwLs1oqbuIeFEJpQCvsPT7+o4dMHQo+5s1J2fmbEa3v5putz/D+6d0wszyl7lne+xRUvSXgo4/E5FwaQZegNey9UHnN81/uLn5tz/ovfFdBrwzmYN+2cbMlp15uP21bK12OOyDFDvAmJ6tC5VaQplZh7riUkQkjwK8AK8QBd9hCC2+/oynF0yk5ZaNrDi6OaP73MuSWo0LXaNoeSTQLwV/91dgi0ioKvxS+mj0VvcYPIXr33iKi9YtYVO1I8jocB0zmp/t2c9twDcZF/odQ/WUZMwge9dezbJFJCQJuZS+xL3VO3fCI48wZfQjHLBKjE27gv+060FOFe9+bvBfHklvk6pebxGJqgod4F5tgf1fWckDM1bnz4SLzYrPbUL66kUwdChs2sTikzoz/Iyr2XxY7ULX8ncWZaAHj6G0KYqIhKpCB3igHursnL1+/15n9acc++Q/YNN6fjvxJAZfMoS3qh9L0WJJJHtsq9dbRKKpwgV4wXpzJTP2h1jjr/v7Noa88xxd1yzmp0NrcW/6XUw7sSO79vk+7yh+0nu4e2yr11tEoqlCBXjRGnMo4Z2y509u/uh1+n48DcPxv2dczn9O68GuKimwr/Dn88I7lFPc/T08DacjRUQkmAq1kMdfjRl8G00V4xxdVy9i4cSb6Pf+y7x9XDs6936aMe2v9oW3h1DKHV57qgCM7N6S1BopGH+daan6t4hEokLNwL3C9YBzjO3ZOn/2e9KmdQxfMIGTN63j86OO4/aud7Ps6BPz3x/opPdQyh2BHlYuHdJJgS0iUVEhAjyvXOFVMKlXI4X0NqkctHUzbugwuqyYz7ZDa3F/+kBePP5sDju4KjWL9GYDEZc79LBSREpD3Ad40bp3USnJSQw+uwE89BAXZGTA/v0wdCi1hw5lRLVqjAhy/UgWAelhpYiUhrgPcK+6N0Bq9YMYa+s59fKb4PvvoUcPGDUKGjUK6dqRLm3Xw0oRKQ1xH+BeZYlWmzcwY8lUWLoUTjoJnn/edx5lKdDGVCJSGuI+wIuWK2rv/JW7Fz/PpavehiOPhIkT4frrISkpwFWiTxtTiUisxX0bYd4+2lX37eGWD17lnQl96frlO2y49mZYvx569y718BYRKQ1xPwNPb12Pem/P4ujxD1Dvt59YfEIaex7O4NyuZ5b10EREYiq+A3zFCujfn3bvvgstW8JrL3J2585lPSoRkVIRtIRiZpPMbKuZrSrwWi0zm29mG3L/rBnbYRaxZYuvNHLKKfDll/DUU/Dpp6DwFpEEEkoN/DnggiKvDQEWOOeaAAtyv4693bvh0UehSROYPBnuvBM2bICbb4bK8f2PCRGRcAUNcOfcu8CvRV7uCkzO/ftkID26wyo2CMjMhBNPhMGDfe2Aq1bB6NFQo0ZMby0iUl5F2oVSxzm3GSD3zyO93mhmfc1smZkt27ZtW2R3u/VW6NaN310Sd17/CI1OvJm06VlkrsiK7HoiIhVAzOsOzrkJwATwnYkZ7uczV2SxYG9japx7My+37sK+Sr6WQB1HJiKJLtIA32JmdZ1zm82sLrA1moPKk7/PyeHN4PBmxb6v48hEJJFFWkKZAfTK/Xsv4I3oDKewQPuc5MnKziEtY6HKKSKScEJpI3wZ+ABoamY/mtmNQAZwrpltAM7N/TrqQt1+Na+cohAXkUQStITinLvC41sxb7r22pbVH5VTRCTRlOu9UPL2OSnIz+Fo+XRggogkknId4OltUoudITmmZ2tSPQ5G0IEJIpJIyv3yRa9tWXVggogkunIf4P7owAQRkTgNcNCBCSIi5boGLiIi3hTgIiJxSgEuIhKnFOAiInFKAS4iEqfMubB3eI38ZmbbgO9K7YbRcwTwc1kPohQl2s8L+pkTRbz+zMc452oXfbFUAzxemdky51zbsh5HaUm0nxf0MyeKivYzq4QiIhKnFOAiInFKAR6aCWU9gFKWaD8v6GdOFBXqZ1YNXEQkTmkGLiISpxTgIiJxSgEeBjO7y8ycmR1R1mOJNTMbZWZrzexzM5tuZjXKekyxYmYXmNk6M/vKzIaU9Xhizczqm9kiM1tjZqvNrF9Zj6k0mFmSma0ws1llPZZoUYCHyMzq4zvA+fuyHkspmQ+0cM61AtYDQ8t4PDFhZknAeKALcAJwhZmdULajirl9wEDnXHPgdODWBPiZAfoBa8p6ENGkAA/dGOBuICGe+jrn3nLO7cv98kPg6LIcTwy1A75yzn3tnNsD/BfoWsZjiinn3Gbn3Ke5f9+BL9Qq9Ob6ZnY0cCHwTFmPJZoU4CEws0uALOfcZ2U9ljJyAzCnrAcRI6nADwW+/pEKHmYFmVlDoA3wURkPJdbG4puAHSjjcURV3J7IE21m9jZwlJ9v3QMMA84r3RHFXqCf2Tn3Ru577sH3T+6XSnNspcj8vJYQ/8oys0OB14H+zrnfy3o8sWJmFwFbnXPLzaxDGQ8nqhTguZxz5/h73cxaAo2Az8wMfKWET82snXPup1IcYtR5/cx5zKwXcBHQ2VXcBQM/AvULfH00sKmMxlJqzCwZX3i/5JybVtbjibE04BIz+xtwEHCYmb3onLu6jMdVYlrIEyYz+xZo65yLxx3NQmZmFwCPA2c757aV9Xhixcwq43tI2xnIAj4BrnTOrS7TgcWQ+WYik4FfnXP9y3g4pSp3Bn6Xc+6iMh5KVKgGLl7GAdWA+Wa20syeLusBxULug9rbgHn4Hua9WpHDO1cacA3QKff/7crc2anEGc3ARUTilGbgIiJxSgEuIhKnFOAiInFKAS4iEqcU4CIicUoBLiISpxTgIiJx6v8BRNxcYGyxzAUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pyplot.scatter(x_data,y_data)\n",
    "pyplot.plot(x_data, w*x_data + b, '-r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It works! That line looks to be fitting the data pretty well. Now we have a \"best fit\" line that represents the data, and that we can use to estimate the value of the dependent variable for any value of the independent variable, even if not present in the data. That is, _to make predictions_.\n",
    "\n",
    "##### Key idea\n",
    "\n",
    "> \"Learning\" means building a model by finding the parameters that best fit the data. We do it by minimizing a loss function (a.k.a. cost function), which involves computing derivatives with respect to the parameters in the model. \n",
    "\n",
    "Here, we used SymPy to help us out with the derivatives, but for more complex models (which may have many parameters), this could be a cumbersome approach. \n",
    "Instead, we will make use of the technique of _automatic differentiation_, which evaluates the derivative of a function written in code. \n",
    "You'll learn more about it in the next lesson, on **logistic regression**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we've learned\n",
    "\n",
    "- Gradient descent can find a minimum of a function.\n",
    "- Linear regression starts by assuming a linear relationship between two variables.\n",
    "- A model includes the assumed relationship in the data and model parameters.\n",
    "- Observational data allows finding the parameters in the model (slope, intercept).\n",
    "- A loss function captures the deviation between the observed and the predicted values of the dependent variable.\n",
    "- We find the parameters by minimizing the loss function via gradient descent.\n",
    "- SymPy computes derivatives with `sympy.diff()` and returns numeric functions with `simpy.lambdify()`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "1. Eric Ma, \"Deep Learning Fundamentals: Forward Model, Differentiable Loss Function & Optimization,\" SciPy 2019 tutorial. [video on YouTube](https://youtu.be/JPBz7-UCqRo) and [archive on GitHub](https://github.com/ericmjl/dl-workshop/releases/tag/scipy2019)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 14,
     "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": 2
}