1b_optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4ohby1fTDHzi"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/khipu-ai/practicals-2019/blob/master/1b_optimization.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3WiD0qVSnGN9"
   },
   "source": [
    "# Practical 1b:  Optimisation for Deep Learning\n",
    "\n",
    "© Deep Learning Indaba. Apache License 2.0.\n",
    "Adapted for Khipu.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this practical, we will take a *deep* dive into an essential part of deep learning, and machine learning in general, **optimisation**. We'll take a look at the tools that allow us to turn a random collection of weights into a state-of-the-art model for any number of applications. More specifically, we'll implement a few standard optimisation algorithms for  finding the minimum of [Rosenbrock's banana function](https://en.wikipedia.org/wiki/Rosenbrock_function) and then we'll try them out on FashionMNIST.\n",
    "\n",
    "## Learning Objectives\n",
    "\n",
    "* Understand **what optimisation algorithms are**, and **how they are used** in the context of deep learning.\n",
    "* Understand gradient descent, stochastic gradient descent, and mini-batch stochastic gradient descent.\n",
    "* Understand the roles of batch size, learning rate and other hyper-parameters.\n",
    "* Implement, using TF2.0, gradient descent and a few variations of it.\n",
    "* Understand the strengths and weaknesses of the various optimisation algorithms covered in this practical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**IMPORTANT: Please fill out the exit ticket form before you leave the practical: https://forms.gle/J4i5wehZPUdggCc29**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "iQ_aoTV805AN"
   },
   "outputs": [],
   "source": [
    "#@title Imports (RUN ME!) { display-mode: \"form\" }\n",
    "\n",
    "!pip install tensorflow-gpu==2.0.0-beta0 > /dev/null 2>&1\n",
    "\n",
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from IPython import display\n",
    "%matplotlib inline\n",
    "\n",
    "display.clear_output()\n",
    "\n",
    "print(\"TensorFlow executing eagerly: {}\".format(tf.executing_eagerly()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SuEFKaqiqe3q"
   },
   "source": [
    "$$ \n",
    "\\newcommand{\\vec}[1]{\\mathbf{#1}}\n",
    "\\newcommand{\\vechat}[1]{\\hat{\\mathbf{#1}}}\n",
    "\\newcommand{\\x}{\\vec{x}}\n",
    "\\newcommand{\\utheta}{θ}\n",
    "\\newcommand{\\th}{\\vec{\\utheta}}\n",
    "\\newcommand{\\y}{\\vec{y}}\n",
    "\\newcommand{\\b}{\\vec{b}}\n",
    "\\newcommand{\\W}{\\textrm{W}}\n",
    "\\newcommand{\\L}{\\mathcal{L}}\n",
    "\\newcommand{\\xhat}{\\vechat{x}}\n",
    "\\newcommand{\\yhat}{\\vechat{y}}\n",
    "\\newcommand{\\bhat}{\\vechat{b}}\n",
    "\\newcommand{\\What}{\\hat{\\W}}\n",
    "\\newcommand{\\partialfrac}[2]{\\frac{\\partial{#1}}{\\partial{#2}}}\n",
    "\\newcommand{\\ipartialfrac}[2]{{\\partial{#1}}/{\\partial{#2}}}\n",
    "\\newcommand{\\dydx}{\\partialfrac{\\y}{\\x}}\n",
    "\\newcommand{\\dld}[1]{\\partialfrac{\\L}{#1}}\n",
    "\\newcommand{\\dldx}{\\dld{\\x}}\n",
    "\\newcommand{\\dldy}{\\dld{\\y}}\n",
    "\\newcommand{\\dldw}{\\dld{W}}\n",
    "\\newcommand{\\idld}[1]{\\ipartialfrac{\\L}{#1}}\n",
    "\\newcommand{\\idldx}{\\idld{\\x}}\n",
    "\\newcommand{\\idldy}{\\idld{\\y}}\n",
    "\\newcommand{\\idydx}{\\ipartialfrac{\\y}{\\x}}\n",
    "\\newcommand{\\red}[1]{\\color{red}{#1}}\n",
    "\\newcommand{\\green}[1]{\\color{green}{#1}}\n",
    "\\newcommand{\\blue}[1]{\\color{blue}{#1}}\n",
    "\\newcommand{\\because}[1]{&& \\triangleright \\textrm{#1}}\n",
    "\\newcommand{\\relu}[1]{\\textrm{relu}({#1})}\n",
    "\\newcommand{\\step}[1]{\\textrm{step}({#1})}\n",
    "\\newcommand{\\gap}{\\hspace{0.5mm}}\n",
    "\\newcommand{\\gapp}{\\hspace{1mm}}\n",
    "\\newcommand{\\ngap}{\\hspace{-0.5mm}}\n",
    "\\newcommand{\\ngapp}{\\hspace{-1mm}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "anRCKj57PRJT"
   },
   "source": [
    "## Rosenbrock's Banana Function 🍌\n",
    "\n",
    "In practice, when evaluating the performance of various optimisation algorithms and hyper-parameters, what we really care about is the performance on a wide range of real-world problems, in our case the minimisation of the loss function of a neural network. However, we are not easily able to visualize what our algorithms are doing because the loss landscape for even simple neural networks trained on just about any real-world dataset will be very high-dimensional (dimension equal to the number of parameters). \n",
    "\n",
    "To solve this problem, we will use Rosenbrock's (Banana) Function as a playground for investigating how these optimisation algorithms work. The banana function is easy to visualize because it is a function that takes a 2D ($x$ and $y$) input and returns a scalar output. It is defined as:\n",
    "\n",
    "\\begin{equation}\n",
    "f(x,y) = (a-x)^2+b\\times(y-x^2)^2\n",
    "\\end{equation}\n",
    "\n",
    "where typical values for $a$ and $b$ are $1$ and $100$, respectively. For this practical we'll use $a = 1$ and $b = 20$. The global minimum of this function is at $(a, a^2)$ or $(1, 1)$ in our case. \n",
    "\n",
    "We can easily define this function using TensorFlow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "IQwjMCKCSggF"
   },
   "outputs": [],
   "source": [
    "def rosenbrock_banana(x, y, a=1., b=20.):\n",
    "  return tf.math.pow(a - x, 2.) + b * tf.math.pow(y - tf.math.pow(x, 2.), 2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "J1tg42OeuaxH"
   },
   "source": [
    "Let's try visualizing the 🍌, first using a contour plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "id": "w18SilLGEehb"
   },
   "outputs": [],
   "source": [
    "#@title Helper functions (RUN ME) (double click to unhide/hide the code)\n",
    "\n",
    "def gen_2d_loss_surface(loss_func,\n",
    "             n_x=100, # number of discretization points along the x-axis\n",
    "             n_y=100, # number of discretization points along the x-axis\n",
    "             min_x=-2., max_x=2., # extreme points in the x-axis\n",
    "             min_y=-0.2, max_y=1.3 # extreme points in the y-axis\n",
    "            ):\n",
    "  \n",
    "  # create a mesh of points at which to evaluate our function\n",
    "  X, Y = np.meshgrid(np.linspace(min_x, max_x, n_x),\n",
    "                     np.linspace(min_y, max_y, n_y))\n",
    "  # evaluate the func at all of the points\n",
    "  Z = loss_func(X, Y).numpy()\n",
    "  \n",
    "  return X, Y, Z\n",
    "  \n",
    "def make_contour_plot(X, Y, Z, levels=None):\n",
    "  if levels == None:\n",
    "    # generate 20 levels on a log scale\n",
    "    levels = np.insert(np.logspace(0, 2.6, 20, True, base=10), 0, 0)\n",
    "    \n",
    "  fig = plt.figure(figsize=(9.84, 3))\n",
    "  ax = fig.gca()\n",
    "  \n",
    "  ax.contour(X, Y, Z, levels, alpha=0.5)\n",
    "  ax.contourf(X, Y, Z, levels, alpha=0.2)\n",
    "  ax.set_xlabel('x')\n",
    "  ax.set_ylabel('y')\n",
    "  \n",
    "  return fig, ax\n",
    "\n",
    "def make_surface_plot(X, Y, Z, elevation=0, azimuth_angle=0, levels=None):\n",
    "  \n",
    "  if levels == None:\n",
    "    # generate 20 levels on a log scale\n",
    "    levels = np.insert(np.logspace(0, 2.6, 20, True, base=10), 0, 0)\n",
    "    \n",
    "  fig = plt.figure(figsize=(10,6))\n",
    "  ax = fig.gca(projection='3d')\n",
    "  ax.view_init(elevation, azimuth_angle)\n",
    "  ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.2)\n",
    "  ax.contour(X, Y, Z, levels, cmap='viridis', alpha=0.5)\n",
    "  \n",
    "  ax.set_xlabel('x')\n",
    "  ax.set_ylabel('y')\n",
    "  \n",
    "  return fig, ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "id": "xFU6ljswmCE2"
   },
   "outputs": [],
   "source": [
    "X, Y, Z = gen_2d_loss_surface(rosenbrock_banana)\n",
    "fig, ax = make_contour_plot(X, Y, Z)\n",
    "\n",
    "# add a marker to show the minimum\n",
    "ax.plot(1, 1, 'r*', ms=30, label='minimum') \n",
    "ax.legend()\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nuOXyVa5COkv"
   },
   "source": [
    "And now with a surface plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "id": "UeMlgRP1TObh"
   },
   "outputs": [],
   "source": [
    "#@title {run: \"auto\"}\n",
    "\n",
    "elevation = 62 #@param {type:\"slider\", min:0, max:360, step:1}\n",
    "azimuth_angle = 102 #@param {type:\"slider\", min:0, max:360, step:1}\n",
    "  \n",
    "X, Y, Z = gen_2d_loss_surface(rosenbrock_banana)\n",
    "fig, ax = make_surface_plot(X, Y, Z, elevation, azimuth_angle)\n",
    "\n",
    "ax.plot([1], [1], 'r*', zs=[0], zdir='z', ms=30, label='minimum')\n",
    "ax.legend()\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Rwx9E5RCssMp"
   },
   "source": [
    "As you can see, this is called the *banana* function because it contains a banana-shaped valley. Within the valley, we have a **global minimum**. Finding the valley is relatively easy, but finding the global minimum is difficult, which makes this a useful function for testing optimisation algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KLFs7xXSlL0D"
   },
   "source": [
    "## Optimisation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JT1dKIEMlL0F"
   },
   "source": [
    "### What is optimisation?\n",
    "\n",
    "Optimisation is the process of comparing a set of items and selecting the best one, based on some metric. So, we can define an optimisation algorithm as a program (or series of instructions) that try to find the best item when compared to other items using a given method of comparison (or metric).\n",
    "\n",
    "### How do we use it in machine learning?\n",
    "\n",
    "Optimisation algorithms are used in machine learning to compare different values for the parameters of a model and to try to find the best ones. For the banana function, the optimisation algorithm tries to find the x and y values that lead to the lowest value of the function (here our metric is: \"lower values are better\"). More practically, optimisation algorithms are used to try to find the best values for our neural network weights and biases. In this case, we use a *loss function* as a metric and try to select the weights and biases that lead to the lowest loss function value. But, we don't need to worry about that yet. First, we need to build up some tools that will let us do this.\n",
    "\n",
    "**Note:** You may be wondering why we need special optimisation algorithms at all. For example, the minimum for the banana function above can be computed analytically. And indeed, for many optimisation problems, including some problems in machine learning such as linear regression, it is easy to compute solutions analytically. So why are we spending a whole practical on this topic? The reason is **computational complexity**. There are a huge number of problems, including optimising deep neural networks, for which we cannot find an analytical solution in a reasonable amount of time. In these cases, we need to find approximate solutions to our optimisation problems, which is what we will be looking at in this practical."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Ag4YyzcSytVG"
   },
   "source": [
    "## Gradients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "4eZ0hZHRysx0"
   },
   "source": [
    "### What are gradient vectors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XK4syJWp0SgP"
   },
   "source": [
    "*Note: feel free to skip this if you already understand what a gradient vector represents.*\n",
    "\n",
    "Let's start with a scalar function $f$, which maps the input $x$ to the output $y$: $$y = \\green{f(x)}$$\n",
    "\n",
    "Here, $f$ could be a polynomial, exponential, or whatever your favorite kind of mathematical function is. \n",
    "\n",
    "Let's now  _approximate_ the function $\\green f$ around a particular point $x_0$ with a straight-line function $\\blue {d\\ngap f}$:\n",
    "\n",
    "$$\\blue{d\\ngap f(x)} = f(x_0) + \\red{f'(x_0)}(x - x_0)$$\n",
    "\n",
    "**Reminder:** the equation for a straight-line function can be given as: $y = m x + c$, where $m$ is the slope of the line and $c$ raises the function. If we want the function to go through a given point $(x_0, y_0)$, then we can change the equation to the following: $y = m (x - x_0) + y_0$. Notice, we just shifted the line along the $x$-axis and made sure the function is raised by $y_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6n0es69z0Lrs"
   },
   "source": [
    "<center><img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArgAAAH0CAMAAADhbK6XAAACuFBMVEUAAAD///+/Hi7DLDvBNUPHOkjLSFXPVmLTZG/Ub3nXcnzbgInejpbinKPmqrDquL3uxsry1Nf24uT68PG7L0LNpavPvL+3P1XQx8qoeZmmgaOhkrcUFBliYmRYWFp3d3lsbG6WlpiBgYPi4uTz8/Tw8PHl5ebKysu1tbarq6yhoaKMjI3a2+vIyuHg4e7j5O+jp82Pk7SkqMuprdCssNHLzeHR0+Tm5/Gdoso7PUyTmL2gpcuLj6uvs9OyttS1ude4vNi7v9q3us/N0OTs7fXBxd3T1ufa3OnS1N7y8/jv8PXHydL4+fvZ2tzW19nT1NbR0tROT1BbXF3q6+zn6Onf4OHc3d7k6u/D1uTS3+kLWI0PdbwGLEcEHS8be74efcAthsQ8jshAjMJCj8RLl8xOlcVaoNFpqNV4sdl5rM+Gud2FsdGVwuGTutWky+Wjw9mz0+nC3O7g7fbe5erR5fISebVQV1vP1trv9voUfa4jhrJBl7oXgacahaCMwc/F4Odhq7sciZhKnKY6mqArk5M0NTXt7u48nZkhkYoklYPI5N8nmXy739QpnXWTzblEqoQ4pXhgt5MsoW4AplEwtnFQwYdgx5KP17K/6NPf8+nv+fQQq1sgsWYvpWdAvHxwzZ1/0qef3b2v48jP7t5JsnW+480xqV/I0Mtlv4Wl2rbY4ds0rVg/sWON0aB0yInD1Mc2sVFAs1pCtVxPuWeFxI+Yx6CiyamhwKY1qkU5tUoynkFFuVVMuVtRvmBewmxfvG1qx3eD0I58woaP1Zmb2aTA58bm9egZTyALIg53zIJefWKo3q+047rN7NHZ8dzg8OLy+vP8/Pz5+fn29vb09PTv7+/q6urf39/V1dXPz8/AwMC/v7+vr6+fn5+Pj4+Hh4d/f39wcHBgYGBQUFBAQEAwMDAgICAQEBAewcC9AAAACXBIWXMAADAkAAAwJAELuPDgAAAgAElEQVR4nO2d+X8c513Hd3aRZGTZli3KHZGUND3SlkCBmqOE0gNSckEOjoGuFhZJLEJoxQY20AIFWlrKUbKKLymybEM4CiXcCCGxQiuBkFbISVzXlR0UH9G/wWvumX1m5pnZfWbmeWY/75+s2dGx9ttffZ/v9zvPk5EAEBCIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKIC4QE4gIhgbhASCAuEBKI271Mjni89bHRKnGNNyBu9zBVmrQLOSmXPd56tVDg3lyIm2Kmi0VZNvSsFPLjxXzFfLeT8qTnO6/mi8Q1zoC4KWZ6vJiX5Wn1DVYLcrkyIptCVuRRnzdelkvENb6AuOmmKOuJ7IScl8ZlS9xC3jcbGPPMIzgB4qabvDymvb+CPCqVZHlKf7clSkidlgvENa6AuKmmLOuZ7LSsmFqe1t9sNe8fcJWQ650C8wDETTUlI8Wdkh2/+yeNQOxJmfOQC3FTjZnimgZrFMycwZMRueL1Eg9A3FQjG5G1KOdt73Pa8ZE7Y5QsOGEgbkqZLpenrRRXysv2yuyE7KzTlstqHqF+jsEk37kCxE0llaJcKMqTWoIwIRtM6O911BFNyyMjxXy+Ik3KxaKV+1Zkmee/GYibRiblfFlRL6+muNVyeUKWS+Vy2agkOFLcKcXnaTk/ka9IY7KV2cpcl3IhbgqZ1PUbN1PcCdmx1LI7OZ1Xs4miklVUZNsrI1wXxCBu+ijLeiZgpbhjjt/70/YSw5iW7haU9VrVngoXuV6dQdz0MSLLWk4wYQpacKzGynaN9RCrBeeqLS5DXBArk2aCMGpUca2ymIpD3NbgbAJxQawUTAdNXctWQUFqFVejJBP9BogL4qQqG5mCFUUnnRUCF3GLZEsC4oI4KctG58Bq8447TXURVybHcyEuiJOyrc07IknVKfVPji5YlcgLzFxiyhoag7ggTsximJbiToy6xFOit2CkuNO2OtkIfRAnQSBu2qgY4k6qKW5xwnbJoGCt1UrySFVd0Gkf2apmRFjmCoibOkY0Syt5Ja6qIXSqNcKOyeP6n6bVXpmRF1fy1n0uiTBPQNzUManG0HJ+XFFyvKglAs7nHazJr4qibLUwllcEr+THrXsmyOUaT0Dc9DEuF0rF/KSi55j6PHqxdUKxapk8Jo+VRsakyfxIyTmBO+oo/XIHxE0h06XShCJmuVRS11p5ebzlTY5afbJyqaTkB9MTpQn7MxJSXqY8lZYsEDf9TMtEfWCSmgfQn0pLFoibbsqlqjRFdsWkEccjaC4UOd9YAeKmmkmlgjtOZApKocE/oJb5XppB3JRTVJoQ+bxLdC36V2mLbp/DExA31YzJxWqJGFiU1C1B/J6FnHD9HJ6AuKmmWswX8u4OVvLeyUKFkkhwAMRNOZWyV1VrMu9VqPVzmhcgbvdS8Vp/lXjfYxTiAlGBuEBIIC4QEogLhATiAiGBuKlmgv9zn9oE4qaYSqH1mZ30AHHTS0ndXNR3JEFcIG5aqRa0TXHJybBUAHFTy6iirWdbV3QgbmqZzstykfPhxPaBuOllwmMuLBVA3BST1lKYBHGBqEBcICQQN1VMc37mOTsgbpqYoJ4tnRogbnqoFGWX/ZlTCsRNDVqHl9i0JqVA3JRQ0Tu8hZTOJrQCcVNBdVzTNrUdXgKImwbKI5q36e3wEkDcFDDZbeEW4qYDZZxGlke7pRKmAnHTwIQs57ukmmAAcVNBYayrwi3EBaICcYGQQFwgJBBXVCa6Lat1AnHFpFJI6+O7AYG4IqJ1eLtm9tYNiCsgeod3pJv/DiCucHTfQI0bEFc0pvSBmu7q8BJAXLGojurhtss6vAQQVyzymrdd1+ElgLhioT6fM9LV9QQNiCsYBVke7/pwC3HFo4xwqwJxgZBAXCAkEBcICcTln3JaDyDpBIjLO0qHt0s2+QgDxOUctcNb6Pa/BRKIyzVGh7er52lcgbg8M6V3eLt8oMYNiMsv00VNW7QcXIC43DKhh1t0eN2AuLwyiXDrB8Tllao6MI4SrgcQl1vK6dul+RxxhfaCJxCXX0oswu1u0wbxarxcmD3r8f1mF4hLFCBu6thtNrcajc16fWXRlXq9vtbYisPiswsLF2yx9ELtAnGLzqXZU8Q1fyBuemiuN9a8bHWlvtlY34nu7S/U5hdmZkxzL9Z8wurF2mnimi8QNwXsNrfW6ktuagZhZXUrkuh7sTYnna+ZUfbs7Dxxi43ztTPENT8gLk9Ux8KuxnbWG+0ra2NlbXuD+OKdMVdbkGqWuPM1rwRX5dzsXKgVGsTliKm8XAz+02ysN+qEgBpL9fpqo9FQV2X2z1A+3mo0Vuuusi+vbu8S36VtztRqZ6SZmpEqXKxRstgLfpkECcTlhunR4OM0O1uby4R4i4srdcXWYPbtNLcaZEa80mCV9C7Uauck6ZLx4Zx/wFVC7myYkAtxeUHv8NIP4226BNql+tpWs63f9TvbrbnG8hoTd+drc7aPztR8M1xJzXLDrM8gLh8EHKjZ2SKkXVltBIyx3uxsr66wdnfWoeop71KYwRmH6TQgLheU9IGakk+83djebMlMl1nWA3bX1+zZx/JWZ4u1szVHzjpLyxQC3mMCcTkgwDm8zbWWdLTeeZwl2djatH2LzXXihuBcrNUuWje3RtOLC6cVSc+eXrDdNB8mV4C4HGCEW4+fZHd71RFqlza3ImwbrNu+2XKjvf8c8zUD/cJC7Xnby5fm5hbma+elC7VT52etDoXzJgoQlwOmfM7h3WjJaiOVVmd91fp+q+1kDJfOnJmp1c6cOWMUFZ635w2XZk+rZYbztUtKh8IMsxdqM8RX8gTi8sCoxy7NG1uOBGFlLa4xmd1t67/LZlvfdMaxNpuxZQHn5tQ/n1eT4JpN3DNmgA4AxOWB6qhLuHVau8S0OxCAjTUzZai3oW6tdt7xkdXQXdDS3edrylpsvjZvpgqX1MpvQCAunzitXY4t1DrYXmlb3TM1R/3LLu6c9sIcmRjUQswrQFwO2XVYu9JhYaoTmmaVoR7uhzjtlND+0Wk1rJ5zhmTiLhoQlzd2tx0lKeazLyHZMBdqoZZpp5z5Kqmks1zmdZc3EDchJkZdv2/TXvnajDmtdcdUdylEcWzemQiQSp53SWjJu7yBuImgdHjJOsJGw9a7qnNhrYqp7nLglkTNOQxGKjlDprhSrXaJuOYFxE0CtcPbOk5jq0Almte6YaobMNW9VHN2wWZb84KzRop70XIV5TDOMTq89sN4d9Zs/ao1vqxVadaNfIF4yYULLSF2vnXY9oKR4s5aRkNcvinJRIfXqjstLq0m/SyuF9t6HrMSoHXXmsGeskoImtOn9Bsu2GYYTtNHHy0gbszo5/DaOry2Sj9PiS3JbkP/KelBd75lqMbm5KySRZzVxxjO2QKu3W46EDdWyHN4bZntcoPDFMHBRj1g0J1tmZc5a2UBs7W5c+dmnp9RGhSXZuxLuBmyQOYNxI2VSWe43bWVEToaIoyNLe23w9KW7zc8W2udUJwzc94Ltfnzc/Pnzs7UTjkz33MuBTJvIG68FG3n8O6sChRsDXb19simX05zkShsLVhpwNmLC+qLZxYuOAbHL9Aep3QAceNlOm/s0myfvxIi2Bps60HXZxG5QDyFc7Y2S9zVwjxR6/UD4sbMpBpubTnCkjDB1sDIdF3ThbMLpxUHiWUW9aGzs2FqChA3ETasvu7KtohvoOGdLswp1a4a+fTY2VnKo5DPh9vKBuLGTtNKbbmt2dLYWfaqLtRqtUuu2epp/yfKzpBB2heIGzPrdXFzBBv6Gm2JyM5naxfPzc24lQfmPTcZVR+LcP0cbyBu5FQKU+a32DZT22UhcwQbW+6J7qW5mu2hBjvnZnzcfH42+HyNCsSNmpI5TrO7ZWrbzsMwvNHUEvXVwD/WuRlyIEznVFhvIW7ElNWBmjG1kGCuyNp6cpY/NrQJi7rLEs2dc6e8kgXPFzyBuFFidninLW2FTm2d7GrLzJUk5isgboSYAzW/bxYSwjxFIABr2psiiwuRA3Ejozqmh9vfMLUVfkVGsJ2UuRA3KoxzeH/9C+nV1lyixW8uxI0Kbd/QZz5nFhKIkmc62NGS97j/U0LcqKgo3n76pRTVvzzYWUnCXIgbGSX52RfSr61SXEjCXIgbFRurn32pG7RNyFyIGw1W3Tb12iZjLsSNgu7SNhFzIS57uk5by9z43i/EZUh1fKRqmwDrHm1Nc+Or50Jcdigd3nFT25Vu0jZ+cyEuK/QO7wsp7pL5o5sb0wgRxGWE0eH9TJdqa5ob06wYxGVCdVTT9pnnAu8Ll0K07m885kJcFkzq4Vbp8KZscDEUmrnBn4noAIjbOcY5vGqHd7V7tVXMDfk0T/tA3M6pWOF2MzVPN7TJdlyNCIjbOTu/I8vyx17sssKtB1sxNSIgbqdsrC6+9KxSTAh+QEKqWY2nnAtxO0Nr7774Im3jzS6iHktpAeJ2hNkn6+JSQitaOXeTuM4WiNsBxoEeadkogRFaUSziYjbEbZMfe2LXPEIJazIn6+pfS7QpP8Rti8u/9NHfXuraqQQqjeinFiBuO/zyz8my/GIKN/hgxWbkCzSIG56Xf/GjSsPhY0huPdEWaFF20JiJm8AuPAnxEz+rd3g/320jtyHQFmgRZlGsxN2OZ7QieV75hY/qHd6/+6enu+Mtt8V2xGkuI3G34xoKSpof/xlN24+98J9//vjjT3TDW24TteayEtmXZyNufENByXLvz+vn8H6m/r+PP/7444+9nPq33DZamrsW1ZdnFHFXu8Hcy/f98+9q4fZP1yXpScXcp4ibgIGW5ka1DGCV46bV3Guv7e/vX1f/eO/f/ocyTiM/81mtBPaIYu5biM8ABlqaG1FNjFlVQTPX7eQrgdm7nVE52Jcu/9U/Km/wOfnjf6GvOF5Rk4XLqXrDbNmMcGiBXR1X25w6kW3Vo2I/Y3LnX/5da5T9sdXJRLJAYXc5upoYwwbEdtrMtXmbyfz3vxJTYJcfVcy9l/g8YNCMribGsnOWMnP3Mg7+Sxmncf4b3KuI+yiSBW/UX8N1z5c7gGnLN7EDASLhtlPczL8tEb/07lPMvS8dbzca1JpYFDP2bGcV0mTu1RZvM/9D/ipBskBjJ6pkgfGQjX6UBRGaBOR6q7i3Xd7DW9RkgbgMTBoRJQusp8MSOsoiAvaDiCs9pZj7JHEZmESULDAfa9QPhI+s1RcbwcS9/Jhi7ivEdWAQUbLAfh5X3+JX+B1dAqUKerLwCHEZmESTLEQwSL5bT0VZ7Fdbxd0nblF5GskChZUonkCL5AkIrf27LHRx4TdHfrpF3KvEPSovP4YxMX92ophZiObRnYZWFhN3a5fdtZee+Umnt7eIm3SeUEIuZsp9WItg1RPRM2daQTeSynMcrC8r4zQ/Zff24Irn91WTBcyUe6PNLLAdcIzqYUm9LCbkEm1XnWpa/Hj+Vyxvv/F7HybuM0CyQKPJ/mmIyJ7y3dCKCyviPQa7pW+Y8IWqdOMN3dvv+u6THyBuNHjlUSQLFDaZ//6N7vF0PW4tCfYk7I6+rdKy9nNf/9FPfvJH/u9bT548efKDxL0ql9WBBYRcXzaWWK/PotxXYU3fD454gV92G/q2SmvGX7I6jPCDirjvc/2p731U8/ZpeOvHFusHZCLdEGR7SbDHIv7gHzRtbRsmaFM036+Y65IsmOH2U8RLwMEK4/VZtDvZ6P1fQSq6G78lf2Kx9dQcTdyPfM/JkycfeKj1M97ymObtU5jJpcF6fRbxFkx6oitEXeyPfk2W5edap8X1ucXvVELug85PuPyUHm7xxGQAVtl6EPneYXrSWOc9Xfiz39MOKvv7lq6JMXD7fYq5H7a/8ik93H4zwm0Qdtmuz6Lf9E4b0eW8jbbxJ89qO30Uqy2vGOJ+RBH3AauY+/LTmraPYow8IA2m/bMYdmvUh26shTp37P7NpzVt8xPEj2Y+4vAdirnvNy4/oYfb+xBuA7PMsrAQyzajW47SKHdsfU4Pt6PT5I9minv52xVz361efEUPt48g3IagyXLv9nj2x9WrC1yWdNeXjXA7RbxmF1d6ryLuex42NlTAMGNoWGaLMW3sbJyXwF3QVc4f+awWbluzWw3b05DvUsx9l/TKI3q4xXMPCRLbjuTrSxxmuhva/6ePy/JImXhRwybuw+9RzP0hvQaGcbBEiW8rfaOky88BjBvGsTl/LY+7h1unuNK3KOJ+Gzq8PBDnGRBG0OXjpGZzLGGxseuyKDNwbJzwfsXcH0C45YBYDy8xgu5S8os07SjTRfr5Iw5xH35AMfeH0eFNnphP3Vlf5mKRZmlbp41ROLeqebci7rcTN4HYifu4qN01wxj/QBfpz6Bo+4eLwc6EbNlj6UGi8wsSIf5zznYSPgBXXZK98Kz8wmKg055axH3oAdcxMRA3SRzQp4/pJnIq446i7UtKy+HZvyRedKN1V7sPu4yJgfhJ5GRJM1+IW911Ndq/oHV4x4mX3SC2Y3yf1fkFyZHQkagbm/Gru7utrgxf+oRfh5eAEPeDLWNiIBESO8u3WY9X3Y01LUF57hnZr8NLQIgrfcAxJgaSIcFDqI3SWBzLtG39v8nn9XA7EizcuoorqZ3f9xJ3gjhJ9PT0bVPdepR94J01o2r7oh5uvTu8BC7iWmNiIDESPvbfUnd5K5qMYWNrZdH6HgU13HoN1LjhIq4xJgYSJGFxjYW+ljEwD7u725vmV1/cXJekiizLpeDh1kPch5EsJE7i4kpS05JreY3hg+wOa5cbWh5dKlSIG31xE1fr/LpvEALigQNxrRU/Q3c3tqxIrgXbNnEVVxsTc9kgBMQFF+Iq0XHFFh07zBl219eWbdaudJQ9u4urjYmh85scnIjrWPurZYat9gLv7vqa7b+AlSK0jbu46PwmDT/iKgu1Vbtzi/VGM1SsbG6tOqRdXFHTjnLIrNaJh7gYE0sYrsRtWVCpEXOzsU6PvTvNRouzSsxWY211VB4JVUZowUtcjIklC2fiSmrcXW5xcHG5vtrYajZbDd5pNrcaq3XidiVL1mP1VF4pgBHfIzhe4vp3fr/ocdIJYAWH4ipGbm0SMloJQL3u4qpN2m0zrZ0e1TplPs+U0fAU13dM7MvDXyauAZbwKa6kyeujpzv1tXX7WmwiH7rDS+Atru+YGMyNGH7FVdhoNjZbc1d3luuN7ZZMYrqoD9SE6fASeIurJQsfIi5r3DN8D3ENsINvcTV2mo1Gvb7kIWy90SCSX4WSHm7DdXgJfMT16/y++o5hmBslIohrstNUaDQaykqt2Wz6lWj17DZsh5fAR1yfMbG33j08fNdbicuAGUKJG4ayFm47/jp+4kofcu/8vnr/8PDwO14lrgN2pFZcaVyWix1UEwx8xdU6v63nSL3truHhu79E3AxYkl5xqwVyl+Y28BXXbUzs6puHh4ffTtwJ2JJecRnhL659TOx+Nad9+/Dw8Ju/SNwHGANxKVDEtXV+337XF6UvKYuytxF3AeZAXAoUcW1jYq/e/U33DA8P349FWRykStxqsD0+QkET1zYm9qVhLMpiI03iTuQZlL9aoYr7kNX5vX/4HoTbmEiPuFqHl0EBzAlVXNuY2Kt3YUAhLlIjrj5QUyRe6BC6uLYxsS/fjYgbEykRt1Jg1OElCCDue61k4c1o88ZEOsQtyTKjDi9BAHFtG4Qg4MZFGsQt6+G2yDzcBhTXb0wMREMKxJ32PIeXBUHE1cfEiMsgOtIQcdUJRhYDNW4EEhcbhMROGsSt5uX8JHGVEcHEdR8TA9GRisXZVNBdmtsgmLhuY2IgSjCrQCGguNggJGYgLoWg4mKDkHiBuBSCiovdxOJFVHEnIluNtRBYXK3z+07iMogEMcWtFOR8ROWvVoKL67tBCGCMkOKWIhmncSe4uP67iQG2CCiuMVDT0QY1gQkhLs6RihHhxK2OR9rhJQgjrvcGIYA1oolbHtEHamJKcUOJi3Ok4kMscc1wG1dNIaS4GBOLDbHEHQt3Di8LQomLzm9siCXudD74seeMCCcuxsTiQrAcd0KWx2IMt+HFxZhYTIi2OBuLpwhmEVJc6Z3o/MYCZhUohBUXY2LxAHEphBYXY2KxAHEphBYXnd9Y4F/cSpzFL5Lw4vqeIwUYwbu4SsthlLgaI22IizGxGOBcXK3DG2vhtoU2xPU/RwowgWtxjQ5vgXglPtoRF53f6OFZ3Kl8/B1egnbExZhY9PArblU/qSzeDi9BW+JiTCxyuBV3Ug+3MXd4CdoTF53fqOFV3Ok8g3N4WdCeuBgTixpuI+5Ep8eeM6JNcTEmFjH85rgFuZB4uO1AXHR+o4VfcSsR7NLcBu2Kiw1CogWzChTaFhdjYpECcSm0L+5D6PxGCMSl0L64GBOLEp7ErY5HcDJkp3QgLsbEIoQjcadG4tqdJgydiIsxsejgRlytwzuSfOG2hU7ERec3OngR1xioiW+nj4B0JC7GxCKDD3G1c3h56PASdCQuzpGKDC7E1c/h5aHDS9CZuOj8RgUH4nIcbjsXF2NiEcGBuKXIzuFlQYfiYkwsInhIFQqRHHvOiE7FRec3GngQtyznOQ23LMTFmFgkcLE4m4xrl+Y26FhcjIlFAmYVKHQurpYs4BwptkBcCgzERec3AiAuhTbFHTxh+wBjYuxJSNypPK9VhFbaFHcgOzBkfYTOL3MSEVdpOSS5O00Y2hY3mztqfoQNQpiThLhah5ffCpiD9sXNZnsHjQ8xJsaaBMTVd6gZIV7gkk7EzWb79VQXY2KsSUDcCZ47vARf+6Y3fXX74mZzeqqrdn6Vei7KYmxIIlUoctzhbWGo3/krPzCGuNlszzH1c96vm4vtR9mQhLjTI6KE26Fe3b6jxEsULHGz2b7jxvrsJFJdVqCO60efKd8Jn7vcsIubzfYPWRH35EnMLTAA4vpw3FKvpy8cPQ5xs7kBy1uMirEA4vpwJMuOr/p6U1w8EMEAiOvDAENxs9mvg7gMgbg+sIy4X/MNZsTFFiEMgLg+nLC8OzQQjj6ntj2DHzC9fcD7G4LAQFw/+g3xckM+d7nhyDJyA5L08PsQcFkCcf0w6ri54z43uWIXt1+1/uEHtXiL1hkTIK4/A1+RzWa/MmwV1y5unyn9u9/14Ic+jBExNkBcCp0N2RgNX8AaiEuhI3FzA2FzYxAQiEuhE3H7w6cYICAQl0L74rYxUwYCA3EptCtuLvRAGQgDxKXQprjHkdxGC8SlwGBfBRABEJcCxOUTiEsB4vIJxKUAcfkE4lKAuHwCcSlAXD6BuBSiEvd43wBxDQQH4lKIStzebBadtQ6AuBQiEvdE6yPvRwPONRyD7ioQl0JE4h7NZvvtH/fnAop7NPzmJKkE4lJgKu5AT7bniPqnQ85MoT/4MxYwVwXiUmApbl/2yFCf6t1QNttne2EgG2LePITkKQbiUmAo7kC2V3n8UqkmHMvaVT2ePUTc7M1QTw8meCAuDYbi9mQPK7uRHVUfH+6xvdAbNMHVOJpFJQ3i0mAn7nHF2WNa+TZnT1SPZg8TN/vSE3oPvvQBcSmwE/dINmskp8ezOdsLoT084qxIdCUQlwI7cfuzWeOPh+2/7I9le4l7/RkKv0FJ6oC4FNiJ22sVEnrs4vVnjxD3UuhDSQziUmAnbtZMZYccv+pz2Tb2yen6XAHiUmAj7uDAQH822zcwoG20YH8izZnvSoOHcj3KY+1HenN9nmF10Pk53QjEpcBG3AF1j/Levr4+Yh12xNmKyA0MHsvlTvT3Hjve45kQDIXf2z9tQFwKzFKFw9barOW67df+UbX1ezjbmxsyehWuYLQM4lJgJm6fR/Ggz6bnUE7tpw2oyvZlvVdtfd5OdwkQlwIzcbMefV27gwNa1nBYzQSOH9K9Pdbf1zdwwvOTuhOIS4GVuCe8fvHbm2j6cqzX3hAe6sseHjzWUgCDuBCXAitxj3mlpS7XHbO6h7WEoc9xG8SFuBRYiTuQzbp3u0hxB+0nWZ7Qo++gI0WGuBCXAitx+xzzYDZ6iSXYgL3YNWCkxll7nwLiQlwKrMTNOcq1NkgHHSmu+bKjxkB+UrcBcSkwEnfIa21GOjhkpLiHj6jC6y/32u/r9exNdAsQlwIjcQezXk/nDNjKZEeUjvAxPcXVirqm8H32FRuZGHcbEJcCI3EHPJu0R61V12HV0n7dygE1YbCJa+UaJ7xWet0DxKXASNxDnmMxJ6xWsNorO57NqcH5uKavq7jhR3hTB8SlwEjcHq+1mfKS8Wv/ULZPGswNHM32njgxoGexNnGtB3wOh33YJ31AXAqMxPVcm9ktHOrPqfsuHOvNZg/pNvdYVQXrK/SGH+FNGxCXAqN53Ky3ase9Krwqh4w1mW11d8L3M7oDiEuBjbhH/FTr8ysRHNVzDPt67HDXF8MgLpWOxR3M5QalQ35J6aB3/qsUxbRyhO1pnaEcdgSBuDQ6FrdHSU5zXsUwlUP+Ibd3SJKO56znK0Nt2JRWIC6FjsXNZXMnjngvzSQ1hPpVt4715Pr7stYzP+E2bEorEJdCx+Iezh4+mqM8lDvol0lI0uDgoBWwh3q7PlGQIC6dzhdnRw95P65r3hN8uYXNGlUgLoWYDi856l0uc3K4B95KEJdOXKfunAj4+z/ofWkH4lLAcVF8AnEpQFw+gbgUIC6fQFwKEJdPIC4FiMsnEJcCxOUTiEsB4vIJxKUAcfkE4lKAuHwCcSlAXD6BuBQgLp9AXAoQl08gLgWIyycQlwLE5ROISwHi8gnEpQBx+QTiUoC4fAJxKUBcPoG4FCAun0BcChCXTyAuBYjLJxCXAsTlE4hLAeLyCcSlAHH5BOJSgLh8Av4NOPwAAAMBSURBVHEpQFw+gbgUIC6fQFwKEJdPIC4FiMsnEJcCxOUTiEsB4vIJxKUAcfkE4lKAuHwCcSlAXD6BuBQgLp9AXAoQl08gLgWIyycQlwLE5ROISwHi8gnEpQBx+QTiUoC4fAJxKTz55JNPvux/C0gAiAuEBOICIYG4QEggLhASiAuEBOICIYG4QEggLhASiAuEBOICIYG4QEggLhASiAuEBOICIYG4rlx5/fbtW3uStHfr9u39K253gGSBuG5cO7i1t3cnc/36wWt7tzI3Xe4ACQNxXbh6cEOSpOuZg4Mb0o1MBn9HHIJ/FBdu31Ku7WUydyRpP5O5Td4BkgbiklzLXJU0cV+TpCs3b2o57tX927df3yNuBskAcUn21YCrxNpr1ms3Mrev791EvssLENeT25kD66XrGdXm1zKve90OYgXienKguarxhh59D7Q0AiQNxPXimpri6lwzSgu3bRdBgkBcLxwp7mtKgUHhpj0Mg+SAuF4YKe6esh573aiJ7RsGg2SBuCR7+0qozeix9eZN1WJTXPyNcQH+GQj21JbDjUxmXzKLuhCXN/DPQLCvxNorbxyola8rd9T6F8TlDfwzENzIHFy7euvOtYODa9L1N7SOg03cA+ITQAJAXJL9NzKZm1eka7cymTduaK/uG+K+jskFPoC4gbhh1XHR9OUCiBuIq5mMNmpzJ3OD+x+2K4C4wbitCXsVKS4nQNxgXFWWatKVO5nrIvy0XQDEDci1O5lbNw/egLecAHEDc21vD3Pk3ABxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAkEBcICcQFQgJxgZBAXCAekiT9P/ngMusBXbenAAAAAElFTkSuQmCC\"\n",
    "height=\"250\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IgxZ4df33XZ3"
   },
   "source": [
    "As the picture above shows, we can use the derivative $f' = \\frac{dy}{dx}$ to construct the best linear approximation to a function $f$ around a specific point $x_0$. Specifically, the derivative gives us the _slope_. If the function above was horizontal at $x_0$, you can see that the gradient $\\red{f'(x_0)}$ would be zero.\n",
    "\n",
    "Notice the _crucial_ fact that the gradient $\\red{f'(x_0)}$ points in the _direction_ in which $\\green{f(x)}$ _increases_, and the magnitude tells us _how quickly_ it increases. Here, the gradient is positive, which tells us the function is increasing to the *right*, and the large magnitude tells us it is increasing relatively quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "l3cO1Vl_lL0P"
   },
   "source": [
    "#### What did we just do?\n",
    "\n",
    "Before we continue, let's make sure we know how this fits into the bigger picture. We just approximated some function (any function) $\\green f$ around a point $x_0$ with a straight line. Doing this gave us some important information: the direction in which the function is increasing $\\red{f'(x_0)}$. In other words, what direction we should move $x_0$ in to find a higher value of our function $\\green f$ - and an idea of how much we should move $x_0$.\n",
    "\n",
    "In optimization, this is useful because we can use our function $\\green f$ to *compare* different items, and we can use the direction $\\red{f'(x_0)}$ to get *better* items. In this case, $\\green f$ is our *metric* and the value $x_0$ is our current *item*. We then follow the direction $\\red{f'(x_0)}$ (near $x_0$) to get a *better* item than $x_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nGPIpQGdlL0R"
   },
   "source": [
    "Below is an interactive version of the graphic. Play around with $x_0$ and build an intuition of how the gradient $\\red{f'(x_0)}$ changes with the slope of the function $\\green{f(x)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "id": "JkZcS1aalL0S"
   },
   "outputs": [],
   "source": [
    "#@title Helper functions (RUN ME) (double click to unhide/hide the code)\n",
    "def f(x):\n",
    "  return -np.cos(x)\n",
    "\n",
    "def tangent_f(x):\n",
    "  return np.sin(x)\n",
    "\n",
    "def df(x, x_0):\n",
    "  return tangent_f(x_0) * (x - x_0) + f(x_0)\n",
    "\n",
    "def perpindicular_unit_f(x_0):\n",
    "  slope_f = tangent_f(x_0)\n",
    "  y_0 = f(x_0)\n",
    "  \n",
    "  x_1 = slope_f / np.sqrt(2) + x_0\n",
    "  y_1 = -x_1 / slope_f + y_0 + x_0 / slope_f\n",
    "  \n",
    "  return [[x_0, x_1], [y_0, y_1]]\n",
    "\n",
    "def interactive_gradient_visual(x_0):\n",
    "  # change the fontsize for better visibility\n",
    "  init_size = plt.rcParams[\"font.size\"] # store initial font size\n",
    "  plt.rcParams.update({'font.size': 22}) # update the size\n",
    "  \n",
    "  plt.figure(figsize=(12, 8))\n",
    "  \n",
    "  x = np.linspace(-np.pi, 2 * np.pi)\n",
    "  f_x = f(x)\n",
    "  \n",
    "  # plot f(x)\n",
    "  plt.plot(x, f_x, label=r\"$f(x)$\", color=\"green\")\n",
    "  \n",
    "  # add a point showing where x_0 falls on f(x)\n",
    "  plt.plot(x_0, f(x_0), marker=\"o\", color=\"black\")\n",
    "  \n",
    "  # plot the tangent line to f(x) at x_0\n",
    "  plt.plot(x, df(x, x_0), linestyle=\"--\", color=\"cornflowerblue\", label=r\"$df(x)$\")\n",
    "  \n",
    "  # plot the normal vector to the tangent\n",
    "  perp_unit_vector = perpindicular_unit_f(x_0)\n",
    "  plt.plot(perp_unit_vector[0], perp_unit_vector[1], color=\"dimgray\")\n",
    "  \n",
    "  # drop a vertical line from x_0\n",
    "  plt.plot([x_0, x_0], [f(x_0), -3.1], color=\"silver\")\n",
    "  \n",
    "  # plot the positive direction of change vector\n",
    "  [[x_0, x_1], [y_0, y_1]] = perp_unit_vector\n",
    "  dx = x_1 - x_0\n",
    "  dy = 0 # y_1 - y_1\n",
    "  arrow = plt.arrow(\n",
    "      x_0, y_1, dx, dy,\n",
    "      color=\"red\", label=r\"$f'(x_0)$\",\n",
    "      lw=3, head_width=np.abs(x_1 - x_0)/10, length_includes_head=True\n",
    "  )\n",
    "  plt.plot([x_0, x_1], [y_1, y_1], color=\"red\", label=r\"$f'(x_0)$\")\n",
    "  \n",
    "  plt.legend(loc=\"upper left\")\n",
    "  plt.xlim(-3.1, 6.2)\n",
    "  plt.ylim(-3.1, 3.1)\n",
    "  plt.xlabel(r\"$x_0$\")\n",
    "  plt.show()\n",
    "  \n",
    "  # reset to initial font size\n",
    "  plt.rcParams.update({'font.size': init_size})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "AucAFvjAlL0Z"
   },
   "outputs": [],
   "source": [
    "#@title Double click to unhide/hide the code {run: \"auto\"}\n",
    "x_0 = 0.81 #@param {type:\"slider\", min:-3.1, max:6.2, step:0.01}\n",
    "\n",
    "interactive_gradient_visual(x_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OjUiG7calL0d"
   },
   "source": [
    "Again, notice how the gradient $\\red{f'(x_0)}$ (the red arrow) always points in the direction where the funciton $\\green{f(x)}$ is increasing. Also notice how the magnitude (length of the arrow) of the gradient $\\red{f'(x_0)}$ increases when the slope increases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "86E1GNVPwIC7"
   },
   "source": [
    "### Moving to higher dimensions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "b1ojWHK2wLi7"
   },
   "source": [
    "In fact, a similar thing happens if we are working with a function whose input is a *vector* rather than a scalar. The same idea applies, but we must replace the multiplication between the *scalar* derivative $\\red{ f'(x_0) }$ and the *scalar* difference $(x - x_0)$, with a dot product between the *vector* gradient $\\red{\\nabla f(\\vec x_0)}$ and the _vector_ difference $(\\vec x - \\vec x_0)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rkqEe17x32eQ"
   },
   "source": [
    "$$\\blue{d\\hspace{-0.3ex}f(\\x)} = f(\\mathbf x_0) + \\red{\\nabla f(\\mathbf x_0)}\\gap\\cdot \\gap (\\mathbf x - \\mathbf x_0)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "r8WVIvRx01vJ"
   },
   "source": [
    "<center><img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAArgAAAH0CAMAAADhbK6XAAADAFBMVEUAAAD////GqKrCLjq+Hi68Jja0NkLIQE/RYm3hmKDtw8f03eGoU2OMi4yihrLm4/CalsavrNJ/fK2SjsWWksaem8pua6t3daSMisBubJRfXp+GhbiSksZeXoD6+v7Y2NlARIh9gbhQU2tLWZNAUIx0grJCR1hrfaxgdqUSExVKbqJEaZpRcJtUeKJihq1ffJtCaI5HbZLy9vpagqZiiq5ehqhLaYJRcIlZeZQ5i8lWhqpLc5FiiqpehKGfxuQ4PUEOdr4efsJGcpJijq6Bttzb6vRWoNFejapijqpDYHI7aYNBcIxGdpJLepVHdI1Of5pUgp1XhqFci6Zikq692+xGepY0TlxikqofJytGepJYiqJekqpCd44qNTlARkdLUVL2+vrd4eHKzMzAwsK4urqvsbF0dXXs7e2nqKiYmZlmj4nu9vRWW1ml4MPH7NoFp1RoyZeE1KyfoaAptGtGv4Dh8+jR1NIeezslikIvi0pka2Y1lE0aMCArlkQ1m00xo0Y2o0o8m091oH02sko2rko1nEY0lEUzjEEzhEA6sko2o0Y6rko6qko6pkk6okovezo6kkdFn1I+qko6nkY+pko6mkY+oko6lkYtajVAiUlGkk8+tko+sko8rEY+rko3l0E6kkI+mkY6jkJLplRSnVl+p4KQuJSkxac+okY+nkZQjVS60rxEtko+jkJQslY+ikI3dztWqVpEd0Y8gz4tUS5Mgk5fmWBurG9KukpUulVxtHF+xn5+wX5rmms+Vz6Knor6/vr2+vbk5eSBgoFgs19rrWl6v3h2uXRtu2mCyn4pPChevVZouGN+xnl+wnmCxn53vnCCynqGxYB1wmxzvGp8xHKCxnqGyn5rwV2GynqGxnqKyn5zwmB9xWyDxnSGtnuGynSKynpIZz9+xmReiVCGvnKOynqMynVUeEaIyWySynpsmVZ1pl6SznaSynaSznKIv2p/smWSzm6Pym2SynKWznaRxXKWznKWzm7Wzs768vL+/v76+vr19fWSkpICAgLnkulhAAAACXBIWXMAADAkAAAwJAELuPDgAAAgAElEQVR4nOy9CVhUV5733zxv///zzsyb7h41rcZ2YqMhiUsAlQIXAoJaUotVlFIjFODrbtzXRIl2Ok6zxEhIQCMx4gZGJUZBEbUgdEeqAAGBsUCRWNUsBmRQo4GkSEcz73PuVufec+6tlTX1O88JRVFUVL5+/dzv+Z1zf+XhLncNwHIL110DstzCddeALLdw3TUgyy1cdw3IcgvXXQOy3MJ114Ast3DdNSDLLVx3DchyC9ddA7LcwnXXgCy3cN01IMstXHcNyHIL110DstzCddeALLdw3TUgyy1cdw3IcgvXXQOy3MJ114Ast3DdNSDLLVx3DchyC9ddA7LcwnXXgCy3cN01IMstXHcNyHIL110DstzCddeALLdw3TUgyy1cdw3IcgvXXQOy3MJ114Ast3DdNSDLLVx3DchyC9ddA7LcwnXXgCy3cN01IMstXHcNyHIL110DstzC/UXWnNk8v+uFryNP9c9yC/eXWHNmzeX5XfN/pZ+VW7iDtF6fPXvutIXU7232rODguZOZ3+icWdN4f9MLgweGct3CHZw1JxgUJdzZwbM9ZlsEOXnarMn8v+nX6e/q3+UW7qCtucHB5G9tYXDw2teDgxmXnS0szWmz5iDP9b9yC3fQ1jTaY+cGT/OYHRxMX4/NsQIDc4L5rtz6U7mFO1hrMiPVWcFzPSbPfp3Gg2nChguEbuUF/aHcwh2s9TqNuHOCg+GMa2Ew/5UZ/QphS+4X5RbuYK3ZNOIyCiZrdrDVqHZWcP+nXLdwB2tNo52VUTBZNqjSBm33ebmFO0hrLYO404JnQb/FhazP8PW6VZro+3ILd9DVwmnBwXPnLCQAYfLcudOCg2dNmzuXWYtgZQaTqaWJyXODg5mXACru938obuEOspozLXjuHI/Zs6YR4pszdxoQ7txpjCqnwcKdDJ5/PXjanFmzJ78ebGGI4P6fK7iFO7hqzixSmMBByd/YwmAW1LIQdy7xmmnBoEFhFvS6af0fct3CHVQ1mb4kW8gIdzZbuLCZvk4u/c4lXhEcbFkxm9v/1yDcwh1UNZvW5UJGoOxrs8mwjKeRvjqLEPvrr0NfcAvXXb1Za4Nplc4ODl5LPprGWk5YCF93kYY7JxiR6ez+vwThFu5gqrlQBkYlWpPZslyIBgYL0Uux2f0/D3MLdzDVLFqDFrlyZIkR7tzgYG6Xo9tx3dWbtTYY6mSkcgHOtdkcVKVzUXud6xauu3qxLFGCxUXnciwW4YK1KOK6UwV39WrNgRB3rofHwrkY8wzmtiowPTjTLIp2C9ddvVkLaeHOIax39lxCqGzhQknXnLmzSUsm8oeF0H6eWcjlWr8rt3AHU1HLZh7TCBedtRCGXaqgXgVirWzyLIolplleNxmx5f5XbuEOpiJjrMlzCeG+PgttxoU7vxYGA33Onkbi8FzoCs2WDrK+LrdwB1NNnhY8e+HCWXMnzw2evXDW62gzLuGm9MNZwbPnzJ41Zzb4nmnQ7nWP1wfArjO3cAdXzZ42F/ybP3n2tLlUsw032LJsOQM9j7PneHi8Pm3utNlrWS9xN5K7q2+Li7jATq1ltHMGACm4hTuIazJ5pALnNzh5ltBxIB4DZOeOW7iDtojOXJy9vm5Fl5NnTbOi7P5QbuEO1poLzq7BnqEgeALTADFct3BdUU+f9cNf1Nzg4IX4ExKEj7LBf0+/K7dwXVDLFgUVL+lvv6g5s+bO5vlHXwgWJk/Df09/K7dwXVErh4/0jF35Yv/6RU2e/TqfPufyr4zNFuaIflNu4bqilnoOGzpspKfnytXXBsYveDZf4PX6NF5J969yC9clVTxkGKjhoz1XvtHPjBdffLbK93y/K7dwXVI/ev1h2NBhYA4B2l06CH5L/bzcwnVNxZGWS9So0V5xP785GH5X/bjcwnVNLYkZNhQew71in1vWH1OywVJu4bqofh41jFNDvIKeczNDT5VbuC6qVTEk47Lm0NExQSvcvtsT5Rauq2rlEK7lkjUyJmjFqsHxW+xP5Rauq+ofi9iUC42RQTE/u7Xr2nIL12X183DEbS01ZHRQ3LJB8hvtF+UWrsvq2SKEcVlz+GivlW7tuqrcwnVdbRmJGC3iuwNlZa3fl1u4rqunQUNRvuUOoF13SuZ8uYXrwnrOquWSK2sjvVYWL306aH7bfVJu4bqwnnmNwLAtbo4Y6RW30r2y5kS5hevKsk65cI0JWrncnZI5WG7hurS8RqBUKzBGeMWMdmvXoXIL16W1ZTTiq9ZqZFD5zwOk/bw/lVu4Lq1VcaOwTCs8h3sFPecOyewrt3BdWw5YLlFDxsT0tz1r/bvcwnVtPfUcIsC0gmPI6KDFWwbTn0WPllu4Lq4Vnoib2l6jRnqtWz2o/jh6rNzCdXE9XWlrloufYN9PsXtlzWq5hevqWuaM5RI1YqTnymL36oRwuYXr8lpsX5aLHWBlrdid8AqUW7gur6VeiIc6VMM941aucCe8POUWrutr3QgefrV/jg6KGyiH4/RyuYXr+nrTacqFa3TQYrd20XILtwcq1gWUC48hXouK3QeMsMst3B6oJUGIbzpbQzxjVv446P6gnCi3cHui1o1yEePCc5TnonXuPWt0uYVrX21ff+rUqc0lWqLKcnNPn0k7k3bmzMdgZOpytdpNWq325mcxL7804g8ut12QksUWuwNeD7dwbaqnL57668k1+VlZKSkpSYkJCQcTwIxPOJAQD8bB+IQDByjtgvHB6TMf54WN9fX1Hes78ZWJr7z88ksuZN5RIOF1X6y5hStUfzt16rPLuW/vSNqxOzGB0OnBY4mJKUkHDpxJS0sDE9QHaR+kpR04EJ8WH58QHx8fn5a2//SZGF9OvTLp5ZdfGop4qCM16gWvoJX/EPh1/xLKLVxMbdhwavPJy/mJiYk7EhMPHj1w/MzxM2doTz2ddjAhMTHt49Mfn/54/2l4JIP5wQd798bvjd+l9vXx9fH1DSGGpSYC/S4YgWNY2+eQoUNH/tJXJ9zCZdXaN9dv1urfBTxw8EBCQsKBNEqt5KC1+vGxxKNp+0+nJpMzY//+5IxkemYkZ+xPPUpYrpQYhHZnQHPixFdeeemlESMcRl3imLJf9sqaW7h0bf/81MnLWQmACQ4cP5B2/GBCwgHIZwG7fnya8tfk0/sTEhM/ID02eT8Yp4n/JqcSMzU5ddc42m2lvj7EpD4P8Q2Z4EtM4L6jRiAUa8MYQn3XcM/Ro3+hB4y4hQuuvt489dmJxMSE4weBYtOOg6QgYTfjtB9bFEsMyltPJ8Xvz9ifuj8jNTU5Y1/q/ox9xEd6xswICfGFJ61d5jmKIohrN/tqKHQy5Egvr1/i0f2/eOFeW/pZflJCwoGDQLXHgWoBzx7YfQBKCT5mqRZ47GnCV9N2p6UmU2Nf6r7UVOK/+1JTM/alpuapQ3yRQbov8rzvxJde+oM9nDtkFPz5SK+gn39p2v0lC/en7Z9vPpSQACv2Y4JpTx9NIJwWvvIiWNbCsMBn92WkJiRkpJLjSOqR1Pcz3k99n/54VO07MSQENwntzmC+NoOYM3wjXraDerknQw4ZE/vzL2p14pcq3J+WnMpOSjx48MDR4weOE0x7hhiEx+5MO3P6A4hn4UkzLOWzyTuSU1NTM1LfZ+b7GUfIx3lqxFmh4ePjg3HekFdesuH8MXAiwxDkuVGjFy3+5Wj3FyncVacuf3ICsAHNtGcA1QK3BUz78W7Ia5msIBXwbDIxSYbdR3lt4l7gr2C+f+T9I0csc2+Mb8iMkBD+CaiX/Dr4QEzwWP6KLcw7BPeaESMXxa74ZRxK9osT7qpTX2aDbJZkWtJr084cZ9Ztz6Tt/uBjyGsJpk1NTt0PsWwyybGUx6bF73s/9UgGOfbtfSf+nfj43Qm7E3anx/ggRaS7kFJ9pVIf+jGUQkx8+SWEa9lz6BD0OWIO7383Z+2J+kUJ96cln+WfOHGQ0mwam2upBOGDBJppMyxeS/As6bEMyzLzyN7d8Xvj4xPjwX/3/iV5396/7Mt4//CRw3vzQkJCJkITCJVSpo8vqWLwvNRXCr5Opr7UeoWPr2/AK4LIO4T/q8M9vVYO9jbIX45wn35+MvvEsYPHDhw8cPQA4bXHKa5lstoPTn+ctpvFtckQ0xJ5gYVlj2RkEM4aH7/36O74jIwjzDh8hBrp6hCEbeEBlAoseZyPj49EKkWI12fiC39E6JYZv0eegcbwMYP86P5fhnCfbgBUe+w4xbUgQyDJluZaIq09vf/0B/Ew1+7fDzEt4baExyYng5aE+Hfik/dmHDny/uH338uI3/se+AjU+h4xiMrYFSLEuPSUSgn/BQW/nqSHiRNf4rHWUcid1Ti+Ozpu5ZZBCw2/AOE+XXLq8iHgtQcPHgdsa0kR0iwrY2Ramxa/n+Ha/RauJZkWjL17Cc3u3bsv472MI5a5Lz7tvfcpn4UqJSaEViBf+fpOkNK8S6mXSBss3Osz8ZWXsGc1EGfmYJ63zCGenqPfGJz954NduE+XnDqZfezgwaPHGbIlFJuGWxc7HZ8M8lraayGuzUh9n9JsRvKR9w8feY+ewGHfBzM+mfTZT+GZERMiD5GFyHlEGwL4lq1pX18f0n+Bbik+9p3hM/GVUajvjrBiuUT9ZrRn3F8H4erEoBbuv19bnZ90jPDaA0fJxDYNw7V0hvBBPM22rKw2NTVtb3xCYsJeS3bAoVni4zt7D2cAPGDme5++92mKGpEqS7YTkKdCQsi+Bov3krw7zucVpK9huI19vsO9Yp5bPshSskEs3Gef7yFUi3LtGW63F8G1yfHJ5JpYMtN3sC8j9TS4AksDbEtmtATLkjwL0SyoPyfTjz49/Ol71Hg1JCQiZAYxQ6iPlin1hT9n+a4UvJpIHnwp5vX1mRDyMifhFbqzGrvAfstlg0m7g1W4z5Z+lvThh1iuhd2WWdFNTt7/JxbXpqYmp+6N352YkEzmCITDWrgW5dlPD8fvO3wYaPlTaOwgLFdGDE75IF7L8l3KjX1J9YZQXTozXoJ7eUfZs7cN7LccPIfjDErhPn3xs08Yr4Xc9uMzrBThY6gPITkj/jSLa5PjE5OSU5msFlx5WdiW8loLzxJKPfzOPkqxUO2SykIiQixTTs0QqeUxflr4d6Kvj88EwLxSH9+JIQEvW3azYdfP+GvoyJjYQXJTbFcLd8N65Cmy1m9Anuqhevb5pRMnjn1x8BhYH2O5LcS2H0NsS8z4DyxUuzd+R0LaPiKvJa7L2GRLc62FZxmHjd/36XucsSOG8Fp6WgwV8VhcSX1pCvYlr9iIqzbfSQvo/oRRKNFaGcNjggbD0f0uFu76jXz6XL+RT9IuradLv/zixMFjx48dxGS2eLYlMtsP/kxxbXJiQvyfk6l1sQy694Di2ve5XHv4MNtj4//yKbd2iUJgzyV41tcX4V2+GSJlUt0QXybrDZn48ghHLJdCY8/Fcct746fRg+WkcNdv3bh1E/SpgDqFvuaq+sdfs8gUwRrbfgyxLZHZxhNe+8Hu3fEf7LOsjbG5lkoR4OL6a8Y7lsdUnY6RhbAH0CL3OXrMQG1XynJnsM5G5ryT/viHYUOHDkdY1qY5ZLTXwD6kwTnhbt647attGxnlbt+4GXmJpTZt3I4858r6aemTT44dO34QZdszMNue4bAtmPvfSc5I3ZuQcDQZrI4xvV6WvPbwe1Rm+x4rp+Xw7Keffpqx9zD7ibOf7pIwfiuTg0HoNsLKlLN42BfOg8HVGpHxKl7+wzBMvmtjjRrtNYCP7ndKuBs2bly7eePGbdSnazduXYu8xlJrt24U+rKTteRU/vkTlNsKsW0ah22J3q8P4uMTj8YngyyB7kXA5bVHkJyWPUjlclkhMYZ0U7qkMpmMxb1Ck6kJvrDrgjUKwne9XxpuhWmFxqjRnqN/HpgpmVPC3bRxm8emjYzjbrICA+sFDdmZerrhs+xEIkfAsy2nI4HcgcP0I+yNT4xPS03OIPRqYVsW13LqU5zbkrV376dnPz1LeC01dklkclkENaNDZOCjFb9F5wwi26W415foMiMzXv+XHfbcYYTvesY9NwChwSnhAs1u2LaN8tHtG7cir2DX1p6Bha+WnfzkxAmUbY/TbJuGJrd0r23qB/EJCbuTk1P3ZVh6v2C2RfJaDNmyNHx2bwalW3oejWHcVhYtpR6QHoyf2BFCsDF9mUZwL7nGNnbGy86d0zBipOfo4oG2suaMcDdsZFmsNcMFlrsJec7pWnUqG6iWn20Rt7Ww7emjiUc/SI0HiHCE8VvAtgzfonntYTSrtfAsMeLPfsoeb0nkJN3KZFL6kS0jQsb13miCj4Ff+0rplMxXNDHE+4+IldpXQz1jilcMpNUJZ4S7mSXctRutIqwNL7Gzftq++uKJD2m3hXttabql+r/YKw4k36Ylkrt0/wyEy8kSILbF5LV4r6Xr8G5axdTHxBjGbyNlETJ62lhsBp5A264vvfLmIxoLVipecuiEBnh4Llq3esBo1xnhbtq4EfpsPXORxl/bXByJXfsMMMKHLLY9QO0w581tyez2g92J8WSv7Tsk3Vr2jeEzWx6VsryWGofjz5L1KfUxT0LwLfDbSBl7Rtg4Lb5Ls26ID8G7ISG+InClFiJ/zfnzIUfGxK1egvwx98dyRrjbWFC7yYZLr82uZIWf3sy/eOIYoVt7c9vTCYlH6Q6wP+8le23hLIHgWzbdok4L6ZjNtGfPvp9MqZaaR2OAx8pYfktOO4ukX6mUZF7ac8GlGqHdQGzfrn1zeEzMQFhZc1S42zdv3rxx49bNmzfTS2UbWYtm2zdt27Ztu4fHhk3btlkEvYHl0U4VkC1Q7XGoJyHNNrZN252QvI/eR5Ycn5F6hNo/9pe97xwFYwdZb721Y8eO3e/sjd/7zjvvZPC6Lcdvifrz4Stnz569QhnulbNvqe3lW/4BfBd0OsgI5RL5rohkB9+QkFdccSrvEK/+f2NhR4W7fuu2rRs3bt26bSv1jz9bkxs2btuwfuvG9Zs3bt6w1WKzaze6KFd4ujwuEbDthwzbwkkCnNsySQKV26btTkyD9zXs3pGSkv5W+ls7srOzsnJ1elA1ZDU1GWoMerJ0WVnZWYd2vAXEnLgbkS6i4L2HWZ67Nw/rtw54Lp1ISKVSGeO5M0TkGhvoI5v0R5Re7R9gZa1f35zVCVTYwPLYDTA3bCDkvH7j1q0bPNZvhOB3K28vgz31dMWl9Jycizk5KRcvIkkCN7dlKAHkCPGJCafJfWTvpOxIf+ut7OysB1qt6XbT7aY7TV833WlqYk2mWunRaqqtBSomFLybw7ZnLXXkz2dZlSfB8C09ac7FPcc3wXpaNKFcwLsiKcW7Pj7SGZOcCnaZGuXVn+NdJ4S7iZURbIaFu43w2A0bwXOb2cJ1fg1i1epLSUkfAro9djDxYlJKEuO2BzC5LbROtjcxIX7//lTgsCnZWXr9rYY7DXcavm74uqEJHZxqRUajVp+b9SHQ72GGbSHt7j3CEu7pRajXghmJec5mTwZhrg/BDL4iS8Tr4+P7ykuuuO/EqCDkz77/lBPCZV+bbYbkuYFU9HoiLtuwadNX0Pc4K9xrf81OAkkC02977ERSIkO3BNti18mS4xNTdqQkfpKdrdffun375te3ydFEjq+bvuZqlSxIr60PWy2Trrayyqysd9/e8c7eTymepeZe+jEYZ/8CKJc1IiItg5dnrY1oabQEeK5MNJ/uayCz3Ul/dMHp5yPfQP74+005IdyNLOFCJOuxiXy4CSXarU7GCqs+u3DxBKHbLw7Se8kOHkhKTIToFrPD4XRiSvqFT7Lz9SW3G27fvN1AjjsNtxu+brjT0IR6Lqrdr6HJHq2tTa3GJ/lZHx7csfsdC9cejof8971PT+/i+Cx34LwWnUhJpWGE0Y5ldZT5+PhMesl50o3tv/GC48Jdu5Hlnhgv3YauAWNeZUf931P3Uki3Pc7e33AwKYkIE3DZ7dHElJScS1n6gts3mXH769u3LW6L8i2kWUq3razJFOnB1Girzc16d/fujAxSrBmpV2i/fe/s2SvpEpJfI7AD4V6hyeLeEKlEJosI8RdZ+n2pjOGVUc5GDEPWOfHT6tlyXLicBd9tqJduRJ9yhnGffX4p5QRNCQePAa89fpSm2wMpx1jZLfDb/aeP7khJ/+SSVnu7+no15bU3abflpVuc3+JHa2sbNaFq/OjtnbvfOXz2ytl9RyhueA9EDFcW8XltpI2sy8e90rBomSzEX8rZ10a6LoZd7ZieK5AfQz8px4W7mQ0CmxCVcpRNlOOO+9PS/HQL21LrZHB2ezHhOCu7PZr4Vs4nube+0ZZU3755u/r2beC0xLSbbnFs28r1W9bIfXv37r0ZeynH/U/ivylqrNcyg+LZSHkEM60QLj2kkmiZ3N+Xta+Ncl1nE4aY/rpDzXHhbmMvJmxCuACDuI4z7ptxOSnEOhnFtkTP7VFWdpuQxPjtAcC0+derb5Vqa65Xo34rwLc49eL5lkMN9MUaPduz3n4rcS/w3P+krtd2IR7LjMjIiMhI5Fkb/JgskShEKhPNh/oa6N0TPj6vONXDMHI08qPoH+W4cLeyexM2I60KNOJu2AQ/h5iwLbXky6QTSSc+JNcbjlvoFkoTzqR9fCAxDbjt0fT0C3qg1QKtvrr6NjlYfHsH9luB7BbhWxv9lhmVuUk73vlPynevpKgRl42MjIyUWyb6dcSVcewbIZNIZVKpCO5pYM5q8PEVPvXRSnn2081pjguXQ7DrIQNeD9aB19Iv2LYJ/iYHhPtsz8WLJ1LOnzh24iDac0v7bRo4TTzpaGJeTnZJ9fXq69W39NoCwmsRvxXgW9Rt8QPDtuDzR6z5qLWt7FFr7c4duw+Tnvsq1mfhifgsbuB8VyKVhYhEzLoa7Lq+Pj4vO8G5Qf2zYcxh4XIJdrtFuORa2XoqddgAAcPajRu/8rCzfjqVc+LEiaQTx1g9t+yuWzK7PZiYnp5/g3TYgsvaapJtbzJse5OmWzv4Fpvd2uS3j8BsNIGPxo927Dhy9sqVoxDlyiNxA/VYvsHlXLDLQiSV0X2/0L41wnVfEvZVgRpe7KhEerQcFu7mjZzVW8vnxFrZ2q2k467dCl2ObUBA2GqtvpQEdJvI9CWQ2e3xA5w9ZQnp6Vn66wX6QuC2JfqaQsprEb8FbHu7ockq37ay+RbNbhG/hSfwWzD11MfH+W/v2Ht213zCQbk+a7fnYrg3RCKTSUVSZmsbm3V9fcYKnLQrPPrnyq/Dwt3EbfSytDVu3rh1+/ZtW7dvBJ0K8O51+9saN+QT6w2JiWTPLeS3aXBvwsGU9CwtUOz1W/qq2yXa0ls027Lo1h6+beXlWxvYFsp2DbT/PmrPTdiljmAxLW7yOCzvoFlXKpJFiqKjQ5iMl826vj6v2HK4I6ZGreuP9w12WLhbuea53vLEZuC2az02bNu4cet69jfZhbif7wGUcOJDAhQs2S17T9nxlJQL2XpCtWAU6PWlBderGb9t4PItD93i/RY3MHzLZVto1pnAR2rUxEgkUhzbOuy5kPdKRDLZ6AhZiJQJeNm0C7IxDMNanyP7Iyw4LFz0Msv6Vkjr2ynhWvXZpYsniE7xxBNEdnsUTRPOJKYkQaqtvl1aq9XeopwWy7dNDvCt/WwLzcpHjx4y44E6QCqVRGPoFhp83iowANdGSqJl/vMB30Zb2JfFur6+fKebC9ZQz37YnOuocDegIa11DrBlkwRTy7Mvngd+e+xEYuJBzJ6y48fPJKbkXMqvLqhmdFui1RZcL9HeYvz2pov5FmXbVjzbMrO1sY7y2/9++N8P/7t8hneA/ziJBPFZpz03UiYLi44WAe8FOyrxrOszEb0/mvUxZDHy0+nzckC420EYuxk1z6/Qp9i1dqvtmcKbe1JOkOu7Hx5LgvaUMWlC2oHEtz7JL6m2uO31qlLtDaDXklJevm3qTb6l2NbQZnHchw/U4tDAwImiMIkUw7f0FPBWoREtkfvPJ5g3OsTSz0DuFaZZ12eiA6br2f96yh0Q7jYQGuBWEjZjnrPn65Za+9ccym2PnTi2A8oSmPtAHk1Jz64uAIPW7Q29llawtobw29ssvxVIb/F+iw6EbXHZLdtvwccnlNsSo0g5VSwOnRoYIBonika8FvZcBzw4OixCJCN7GaSsdgYW6xJrafZxblC/20HpmHDXb0fWycivCO0+x38Ppv72Zg5FCceOfXg88ZglSzhIuu3xpLdy9LDXVl+vqtHeqGbYtrTKGt/ylEv5lh51jRDldqhDxeQQTxSJZsgQwo0kC32eGgKOGxkhksp9abqNhnNeuIfB18fuUxj6Hyw4INz1G7et34o9JUz44svWo8N+zE45b/HbE0l0bwLdmZCYnpdfUF0A++2ty9pC6Aqtulp7HfHbHuBbVnb7CO+3wG0NrbTfPnr4qMg/lPBcMFUzwsb5z4e8lb0aLDQRtyWHZL5oPp3thkRD/bss1vX1tfHuEcwY3d/axBy5ONuwadtmvAY3bOW/Ptu01abtZn/bk3ORzMDITrBEmm+pLCEx71DpzQKW21bX6KuY/IscN/S8fGvj3jJUqXbyLeSyxlrosw61OBQainESkWh+JKYw7GudhWUSmT+peuKcMlb/Ltyv6/OKnTsk+lubmIsPdt7AmxvYqNvPcy7SuiXWyhISqd4Equs2PbuAcVvKb0v0JbCKSf2WlPYu3/L5Lem0xkbSbYlRpBJbPFcsDp0ikcz3EfkjGhTyW37+lUpF85lsF5BuJId1GdP9I4Zl+efwOORn1afl6qP01/JludvxHs150Zc5J85TaQLZm5AEZQmJZNcXe5ToS+BVMia/ra3C8u3XvZTfskejEfr8SQzsuJrQeWJVhEQa6S8SzeDYLsK3fANWuyRaBHX4skmXxbq+9uULnv3rZmn96uYlq3OSSN1Se3iPHUxMYLLbpLxsfXVBNYGtcW8AACAASURBVNtvC7XagmrUb4F69X3Ht61sv3343w9rTZByiyZCfksOP/k4aUSkVCTyj3TAc8lJe656vghy4OgQ9p4JmHUBL9g8hvavPb/9R7h/+5xacSDTW7LzNoHy22MpF7JLCq4XXK+G+faGtrQK6qVh9ydU1aB8y5vftrqSb7l+++jho2/rLI+/fxJD6XUeNMQKqSQyMlLq7ztfCnuuXZNQu0ziL4Wpl8x0IzE9DCG+E+3oGutfK7/9RrjXTqZcPM/x24OJCUR2m5iSk1/CdtvqgusFpfoSDjfchvPb0qoe41ur2S2zUkZnCSYTpOPm+VPF4nmEcuGp8JGAZFfq7+8vinbEcynflUr82emuVMY+HxJKdX1fsb1/Iag/wUI/Ee7fPs+n6RY6M+FYIvhvYnpSfgHhtqxRoi29DvvsTaQ/QXsdzW9dz7fc3gTEbYnRZgBeS43HMSTbIlMhlRCbIqLni+b7R9vHuZYcOCJaNJ99lRci43Tv0v264P49Np/AMGpxP2oT6x/C3X4yJ+m8hRPofWU7Pzx4ICn9gr6A67bVBTe0cHKL8duG2zertD3UnwDzbSveb2G3JYYJvj5rDggV44dCSrhuBNCuv2h+NJI02OLDUomIkzWERHPP5LWwru8rtp42NrofHRDSL4T7+QVACV+QfnvsIHMvh4SDiekXSq4XUMPCtwUl+gJIqwjfUvltbVWv9yegXksPA8G35PxmF4tvWUMplchoxJX6+8+PRhNeqwwsQtI1WTQ4O0TiI5FQUyqNplg3JMT394i74qsftYn1A+Fu/zL94kWGE45Z9jkkJqXnlCJeW11QTVGCsN8Cry271cv5LZ/fPnr4qA1ehcgLwDAuPadIpBanlfqLiLU1oYF6bpiI3TcmlUgkkpAQ5VTx1NCp4tfEflNCQiSSceA5oiaOwDAtOoeM7jd3iuh74Z7KufjF+RMnviB7E6B9ZYl52ZDbWkaVvvQWrVVeviXz2xtlDvHtQ3v81jrfUnRrZBj3+0cPFyF8C88oiQT21Rk8a2v8iW+kej6TMkRLJD5KkktCQ2fOE4fOnDdTPFUMpp9SKpFIAev6vmQT6Y7+Gfn59VH1tXDX5udcPGHxW8u+sqS8SyWo1xZUF5QS67t4v0X6b7VVfZDf4vwWXJ/BlHtehGdceijVUpbHzvcRiQQ8F/Hd6DAZ+XGcRK6C1+mmiokx02/ma+T0k0vU0pCQCRNtOn/BC/kJ9lH1sXA/v5BDZ2Af0mcrEXlC3icl1QXXCxC+rbpcytEqD9+S+e2d0n7BtyTbGhjGffT9oxiUbtlDIZFybNVfJJIjXsvHvBLwWolEruL+nZg6byaY4qnQkAPftWkncL9Z+e1T4W7fA+j2Iiu7JdbL0j8pKShAvLa6oFqv5aZiAn4LEtySGizh8tECd7iIb8nxqAz6rEiEsC2S60qiua463wbipfJciUQiDSfeR8zujSA5d2ro1Jmvif3o6TdRMs53og2Zbn85IKQvhft5Nshuz39BZbc03x5Lv6AlVxi4fFsCrTjcZjstyrdkgttgMDl+Ppjr+JYc7W3M467vYzBsy5mq+RI0U5g/H2hXsADlRoeJ/Czvh7AIwbpimnVniqfOnKqQSGww3d/0k5XfvhPu9pM5oO/2PMdvL6Zf0BcW4Py2RFtaLeS3yPkJRI9CXY0VvrXjfDDH+Rah3O8f3QtD1ISOKTLUdSMi58/3hdbWcEMiCYny92YRAtt3KdJ9zcK6U/1mTvQZZz3THd4/jh7tM+F+np0DrsnOU2tlFN8eS7+YX02olsu3BciKAy/fss8Hq2kQ5Fv+88FcyrfkrG0DH7vIkacU+4n9/KaoAgOn+E0BHy0zMEDprVSCOUMikUZzaRcYr/98//ky5GmioiUyv3lisUjMomZrrDtz6sypSh8fq6br+SPyw+yD6iPhPjt14SKlW2qnA6HblPT8kgKs35ZoSwsF/RY5P4FOcE1aYcJFuJYegucnOOC3jx4+6jI8bPvBZCTrfhjQq5+YHCjjqlSBUeF+gYFSSZhUIpKAKZFGWwg3Mhp05KDEG02zrTKA+5441rVwLvDd1yb6vGAl0x0V2x/C3L4R7rU9OYnnub0JB1NysgsL6Q4aFt8WlCLdi0iGi/ItlSfUVFkjXBzftrqYb7saG43txspao7Hd1NbYCXw3LwrDtZgZJSNQVzYf3PEBDJEkmqLcaOl8//ls4pVKFPOo9wgI5L4f23WJTFds4VzCdf1mjLVy5E2/CHP7RLincnLOXzxvoVvCbxPzLhQWFGL9tvRyCbyfF+O3eL6l8tsya35r6/m3fOcnCPtt+w9Go6G2xmgEnxkhdrinRpiWZ6gkEnbXAnDfMIlECvoXo/0lUNIgCbG8R7i/WKFWx6jDwvlYlxx+bNZ9bdwLwpluTD+4NUQfCPcfl3NyTsB+C1YdTqQ3lxQUFGL5tlSPaPY2mivg+JbcY1Zj6Knzb63wbZvJZLxrAIplONfQ3kUz7vdtRUokveUbaKoLShIdJpFER8si54v8/Um6VcHvIFr0K6L+Rx3Ow7ozSdZlk+5r8nGC+4BH9IPrs94X7obzrLUy8tTb9PTLBTx+W/VlYYH9fsvqUNAJ+i1+CPCtDX7b1mg03DUY21ByMFro4ft7aiS75Z9yiQyXH0RERkdLJGqJOsjLUyQJYX2/+tVf0fWqvzDrwpnuzNem+o0dOxTHt9QcvQX5sfZ29bpwT+WQHTVfQHx7Ii+/8Cbht9Ucv71e8HctqlZ+vkUJF/CtqQbBW5gIePhWwG8F+bbNaDAY2tsgvoXS26729u8tnluEZVqeqcKaLl3qMIlaHRsW97KY+f6wX0H1ajgv67L7F4hMd+Zrr4wR2pHW97DQy8IlO8EuQuntiYMn0i+UFPL4LVhycMBvuT24ZVi/dbRHAZspEG7bajKUoT7LysXajFA+1mEz5RJjkoSvPzdCFBQplfgFhnkGLfb09iMJ91VYuL9axEnDUNZlZbrTp4z5PUq39BjS5yu/vSvczy/lnKdUy/Bteo6WYFsM31ZrSwtu8usV5dubX+PPUDAZsIzrkh4F5hrMaKg1trVxdApP0mWNjRbH7crjY1r8mCmTIFZLlleYREZ+T0Cc12LPAO954phfsUvEeU9WvoBmuq9NnTiGv2Osz08T61XhfpaTcp7Dt0l5+YWFPH5beJnpD0PUi/Vb/j1mehPGcbFs69AeszaDwWAyIXkCNhnrgim3Q83LtJgZqgmdjltJi4gMC5IoLK8P9F8c6yX6H45wY7jviWa6oaxMd+qkMfxnL3j18WlivSjc7V/mMGsO5Llgx5LSL5UUVuP5tkB7GdUrhm+vW3oUMHxLJbgNeh7N2t+Dy+HbNlNtnRGhA4RuH1lc1tQIU66Sl2mhOU8Tygy5FLMnQq2ePo/9PWvUHN3+6lX0fYX6F4DrThnLuzdi+MqfkJ9wb1bvCffzSzlJ51l8eyI9R1/A57dEZ4INfmvjGWH6Wzx+i/JtkxDfsnsUGuuMT4zf8Ca4+HXfNijL7XoQg5AsOkI18FBKJVzHlQRFga+w/v0XcYX7K+R9OayL618Y+8JQlHCJ0cc3newt4f7jsxzSb79gTk64mJ5dWl3Iw7el2pLr1fbxreAZCg0GV/coANXWNbahGhXgW4py26DP0oWzXLZmySGWc5rGooPioK/TGS4iXPTduakumum+NnEsT7owqm/D3F4S7oaTZApm4dvE9JzLhWy/tThugR7aa4aol99vkQTXMsq4kQLisnw9Cji+bTTWGusacQkudrCpocsA0cMToSwX0Sw1wiXw/ghpWGwgV9sYx30VfX9iCmS6M/1m+k0cOwrHuH8YNvI55Mfci9U7wqU6bxm/PXYsJV1fUHgT8G0hyreF2N05NuwxQxJcS7XVoprF8W2r9R6FRpOhthGf4fIQbhd7tLdBn+cFouxJTMGhsJhudMi6MNxrucJd5D0T+X+QUzjTnTKGp3chqC/3/PaKcD+jehOY7tsTeRcIr2X8FibcwtLLnBPCEP3y+a3APXobyqoQv3WkR8FUV1bXxDn/1orfIuzQZiD49vsuMJ/wUC6qRNaYLqFIVyILDFKhXw/VcOMwtWK0f0DgTPT/FCoWznSnjMHfI21ULPKT7r3qBeH+48sctt9+mJJXSmm2mtFsAe23JZdLbnLOCEP0iiVca/cwq2E0y8+3VhjXVFvW3oYkuEIa5fItRbkPoc+KAlDyRFSIGUoi05XKNWFhmK9qQpUcUgDvG+jv7++tQv9/sOuime5r4/DhwsiVyA97EAl3+yfNF8/DOVhS3qVCqjB8W3q5UNhvIQ0jfMtLuABwtSaEFTCDP8M1GSpNDvTgYsjhUVsdpOLHSJaLaBA/VJLoCIlMowxS4r++CNbt/8yn3l8JtBuIsi43050JZ7rjfo8y7tBhf/Dqu/v89rhwN1BpAuO3Kela2m1vEqplEe7lUkSlDvAtQrigyI4FmAnwfIv12zZDvalNuEeBT6WI33Z932WAP8sLYDMnokDeIZPINKFhnviv+wdBsPA/csv7qxQzAvyt9OoC1oV7F8Zi28VG9N2tIXpauJfTcy5SfEtWXnYpx28hvi35sop7Ai5mYNbMrPAtqVzDLQd7FL4xGIwN+DUzWwgXYQeQ5XbRjPuw6zGLcu3QbahGrAkNDApEnwcjNmwC0x72qpxNtYFKeYB3QECgrZmu38xx2M4Fzz5rE+tZ4W7X51xk+S1ht7x8q9VWV9vMt5weXGyPArsP94bBoR4Fo8HYiDgwr+NyNYq6LTEM0GNzkb99dEtqViEDm3Q0nuvQr4ER4DV+/Pjx6kWvLlqkFnn5azhMKw5XBgQS2sWTLpLp4hfRYvoKFnpUuJ9np19k8a2FbjF8W3i5tJrlt7gEV8hvBQm3yUK5GLKFM1yW35rKDN/w9ygwhIuoV4AcCC03Mt2N5kfmx3n28K0cnF4nlckUynCNSirnI9zF48aPHz/BKywiYnqoxt8rUINkuFHeAUqV0tvfW6B/wRrnjuorWOhB4X516hLbb9MvXLa4LcK3YK0McVqkQ8HmHlxcmWpQvhXqUWg01Bnb8H0KPAOnUfxoYfzW3PX9BZG17Hb6dCVxrq1MLlMqoyy9C+J1auS1xFgTCwxXFBQYIQfKVseBd5/H7VUIDAgMnBfo7e0dOJMn07WcvYD33NF9dEBIzwl31cmclIsQ3yamXzoJ+20hi2/B7nMH/VawR4HdE1ZrR4bbVlNnaOTpUbA9w8Xw7SOCbY3tBN+aH5m7zI8688Q82a1KoZTJZshkshkRwGHRrwfND8GmuF6eE8ZPGB8rmk44bmhgkHcouaIWyqXdgEAV+K+3twqb6UK9C1jO7aNkoceEu+rLnJzzkN/mpOdDjMDl2wJtacH1AjzfYhmXvwcXm+FSVVdna4ZrMhhMKNdaJVxUo3zjoYHxW7O5634Ypy9hulIpp45gViqjMIqmhihocVhYFPp8QNCE8ePHj4tVRZGOGyrypN4fSXDFSqVSJZ6nAsarEsx0sZ7bR7eG6CnhvnmI6b0Ffnsx/UJhIeK3DN+W6AuRkxRQ9fKvmQnyLdykUCZAuJYehSZjWV1bq1AfrpOEC4ax/XvKb81d5id5KpIUwpVK4jRbqUyuQtSIG0r/xbFBi8MCuISrHj9h/HivME1UZEQ48YzXSJou0AxXNUMZDj6KAwICAv3QTFeQc4OWIT/9XqgeEu6Gizm0aoHfJjF2S/FtNYtvSy5Xo/d4wBEupgfXdsIFmq2rs77PrNFgwKcI9me4PG5LKZfx2+/N3zZLFHJZtEQSIldg/+XHjAUaeqwRqWODFsf5/wfzqoCgiQThzgxVyUjH1QR4Qp04aL8CyBeIxwEBAS/7cTNdpncBl+d69sVNJ3tGuJ/lgN5bhm/T0y8X8vBtAS5NsNNvbclwLZZrhW/rDAYmWuDbZ2ZHlwLKDkx+a+h6ZO5seVCZlfnBzp2HYmQyhVypVIarposxKrUyZLJA/zCvcq8wf1Kdi8MA4XqNDtWER8imk9+r9rcQB8q64lBvhZK671q4t793IP7shdfGjEIod3RftIn1hHC/+jIn5yLEt3n5kGbZfFtdcL1QX4Xc44FPryjh3saco/A1lm/JHoXaRl7Gfdja2lRX1oh6LToQfeI1yj8e1tbq3k7avfPt3NrabwnPvUd1eImjVAqFXKFQRqnwZKuxOC00oiSqBZp5AXGLg2Lj/APDYwHhTohVaUKjGOaI8lLC74Ow7jyxSq5kHgcGBHhTOmZxrh9m969nH8BCDwh31b0cIr2l+DYp73Ihlm9J9ZZeRtwWnyjw+60g4XKrDGFbZp/ZHYOhgb1mhuVbZzLcNkOlLmHn7p07dZUtnUYL4z4y/xDD1qdYNV0hl8sV4SrbWFclVVCUELY4KAgQ7gQ1+LsQHiGj2UMUxvJyDOvO8/MPEDN9un4BASP9LKxLca7fmKFczh0S1/uw4Hrhfn6p+SLEt29dgDVL8y1DuJdLkVNwrfCt3T0K7HXemkaeDLfOUNeGeiuOcXkGolLIb9t0uR/t3PlOVm5lreFbOr9tb2cY19z1/T18j5dYFaWcThgwTLWYIZbJmGeVgHB9g9ZoFgDHpd8rbHEA991R1hXLA1SWPl2xtzfZCQllulPGIJbbB7DgcuF+folYLaP4Fmzi/TuWbwn1Fn6Jvc8Dol/eDJfnHAU+x21tam0wYBm3rswksNcM57c2ZbjtuVkf7dz5p6ysSj3cg/uw6/tH5q6u2m8hz+2KQVVrGarp4UqlXKHEJF/MUEjpR9MJwiX+JigiZNSzgUEif657o6QrDlUpA+jz9wmvneLt7+0H9y5MGsOl3D4Ic10t3M8uWlZ5z3+Rnl5aKMC3JZcLWWfg2rDHDOnBtaFHAWbcptYyDN/WGExojwLvQBSKMm5bbua7b+/8M3BYI550gc+2t38PeW4HzxoYOQgnFauipisUiulRGMfVLNBMl1AEDBbNxnv5E89FyDXkO8Qt1gR4I++Lku48scpboWL1LgQGBnh7W3oXJiJxbu9vQHOtcNeezIH49mL6pb8XCvBtqRY5l5E3URDyW9sJF3hrYw2Hbk11tQ2CZ+FiMwWc39bW6nKzMrOyMjOzKiuNguxA+qzRwrhd5rY822g2VBw1XRmuJOmB48tS4vprCiDc6aJY8FgppxxXFRQQGuWP6STDsW6gQh7I7tP1A9r1o3oXxv6e27PQ6zeddKlwwWrZCYZvL1rCW8pvWXxbqC29Kci3VnsUuOco2EC4wGfLmmC+NRnqTLg+BR6+xfltW8uTyszcTKBZwLCoSrED+KzxWyjLNXfgKZcnSYhSKaPAYLmvRqrQLNAoxo8fHxG6INZfs0ATFSEnvxa2GFy5IZTLx7ozlQGBnLMXZhJNDYAaZiKnOY7obVhwpXCX3CP9luTb9PTLJwX4tgSX3trlt9Z7FFDCBVxb942Fb28Yahr4+3BxGS7Uh9v25IkuMwuMrKzKB0ZE0Wg+ZmFc0mfbas2w5y6yde2BGSqVUqFUKiDvnSEXE4Q7CaxCqEI1ikjKcWPnA84N8EbeA0+64lA/hnWhdTS/SQEB3lNmThnHvVfE8F7exuNC4W64CPyW7mJ86xKXEVhnKJRcLmTf4wHNcDF862iPgmWfGUm5lFK/KTXw9IUJjEcP22r1lVmZf9q586MsXTufRq0Mymdr4VzB/ARLuYjXokMZrlAqVSrysZQgXODSi8NIxgXP+6sJRw4IwKfDCOcSXQwBShXauyAO9B7pHfgCN8716t0w13XC/TyH2MtLksJFIk3A8S3pt/j0ljdRcE2GSyu30kS7rSXDRf0W5du22jJdZlZ80s6PsnIbUcK1luHChEt7rImV5XbtwvR/2TyUihBluEqsmTQOEC54JiA2SqOQRxDvEOtPki7m+gwMDOcSuW4AketyehdeE/t5j35hFCfPjfkK0UQPlsuEe4rsTiD9NiX98mcYvmXOCNOWQqcooDpFCNdFGS7Vh3unrrW19Y6hBO0Jw4824xN9VuZuwmFNqEaxOrXJb783mw1wrvD9E5RyMf4qMMRKhTJiQvCE8WKSa8M00yMA9S7w96L7GgJ49gZxe3WpqQwInAnvR2POXRj7wpAhQ+DjHIf0apuYi4T79LOc5osM3+ZdAGzLx7eFBfj01gG/FexRwBEuybZlTTfK2NkClnHbnpRl7nxrR9LbuU8qjbU1BmNtfcXVc1fPnvvo6rlzV69ePXfu3JWrmVev0J+c+5R4dJX47GpF/V2j0dhiNLZ3UeeEwRku7bHt35hhz90llNPaNvwmTBg3bgL52sCgQEWknCBcEfXdKqU/3/sglEsOZYB3KO7chbEjhg4dMuT3IyyduUsRXfRcuUa4RAzG8G1efiG3qqEMt1BfCPstlm95GJf/HAXbMlxKoXW5dRydsjy3rVZ/bufOt956a/fuvwA9nruaWVlhqDUajSZ7zgrrbG8xglFfAaRccbXi6pX6FmNL+7ewx3ZVwv7bdZdDuULuyjMIwlVIyK/GLVYQjusfy7xaHhi4QIOqlp9054kDZ6hg0gW9C0TXAuG3o4YMGUUa79DevOmkS4S7PT9v1668cpJvc/Iu/12Ib0sv49NbvkTBol/cOQr2E25rXVmdAelTAOo1ZmW9vTNpx1s7MnNzK8uMtUYjT4Zr9SwFhBwodmgDSq6vr68/d/VqfX19i/mR2dgF5wpdeTx5lc1jOpHhhkZJiC7cqCCJXGEhXDDEfNdnoTyZLrHyG6Cah/ToThlD8e2IUb8fMmrY0GEjezHMdYVwlxaVX4jLj7uUtyv9xPmc5tLCzziqvQntMdNqC1zAtwjhfm0b4ba23imra3pY8w1st4asrENJOxOysvT6WgPSryCc4eJ1ilItZ83M3PVte7vRUF9fX5FZUX/37nfEc2az+YnaOb/VTAeGKwaJbhSZ3gLHDfCCXhEQGEhmvqhuiYFyLjiLQenP3Y8mnvrKCwzgDh0xZMgfe/M0MRcI9+eiuJN7Tu45+deTxet25a0rFeRb7Ul+v8UTLjbDtd6jgCPcr1vvGGqbmlq/JlfPjB8Bh333oyxdZW0bh28F+xSs+C2eHbgZLj0NXe3AhSuuXq2ob+ky58kRJdkzxCDDHU98RwjxTkFhCk3oYhH0HlGBIwXyYpRxqREoV4q5PbpjWWnuiCEjRyLq6KlyWrjP1l06efKvJ8n6snnlhZMCfFtyuQRzHzPbehRckOE2Guoa2rS1uqysTw5lHQIOa3Dsfg/IsOksBTRTIHrEjNTj7vb2u3frr+RNV4ZHqRz02wWTCMMlHyuixQs0oiCFJiDoP+DXyFXezEobqlxe0p2nVAZCrgt6dP3GcdLc3rvppLPCXVX0JfDaPSdPntxzMm7xjx4rmr9ELs1ozyV0y+u4iHp5/Za1zwxxWozfNjW1NtQ+yD106MMPP3z7w6zcSn0tbq8Zpk8B47c29YXxswPbb81d5kqab8nRLFKFhyvCleHKqChEU9YGQbiT6FcpQsI1Cq8wjVrE+j6l0pt7ni4nXcBxLsjB/Kaw9wDPnDR2KPs+v569BQtOCvfZ4vyTTMXFgfXqNxdfhhUL7TGjzrOzi3CxPQr2EG5jbZkuKyvrww/fSvkkM7eyxkR6LLM/x6rf2s+4qMfy+y3oEWujH5vBfEBRrgp0IiiV01Vi1Fd5B0241IiSKhWiWP9YFfv1a1TeltcguuXPdDWh4ikqNudyG8VGruylO1Q7J9ynsZTfAlaIpbqJrwXtwfptqR65r44dPQr2ZrimWr3uUNaOpJS3P0nQZeaeO2Ri9eFSuyZ5M1xsn4LNOx+Qic1w6WmC/NbcZW6WsP7dDlcoFOHhPPt4OAMiXGrIQiJjg7jnPoeHe2O6xGwi3VA/bz/W+WLjOGfn9tZpYs4J9zmL3+4JYu58tapcy3Zcgm9BnIDjW0Sxt3GMi7nXAz7DvaOvzM06lJiz81CWrrGmXnfuXGVdU1WZgbPPrMHAyNTZe/YiOkVdVsBvzV3fGx5a/NZs7iI6FjieGqVUKKLCVda4FyZccogVEfOD5nFfp9B4w88huuXp0yVH+BQxxLmhnD1oQ3vpTCanhPtmEMm24LIsCLrv1ZJyJFMovFwq6Ld4wuX3WyTBrSvT52Zn70y6kH1IZ2pqulNTpsvMrawjaKGGk9uCYWjj6VPg7cO1OVNACBfSMuq35kftj80sz23GnoyvCVUpFXJluMDpIGIW4ZLfI5NHcM9c0IQGhqtwXWKI62I4F/TeTPKDOHfSGHZv7vDeuemkM8JdG8T4bX4QC21WlVOaZfaYXa4q4PAt4rQYwsVmuF+zCNek13+0Mycp50NdJX2XhzpD7rlcYx2dKRjK2jBnKdwx2OC3NvAtRqeoz1rx3BbIb81m8+MiAWfVhCsViunYDl2ScP1Yz6nkEUqJivvKecpQ75mc98UoF2VcavpNmWnhXC7m9s6eX2eE+1wclSecbObu3FiezvLbkvwSzn117PNb7prZ7YY7dXp9ljon5cJOnb4Wcl6DLvPjyhooDgObyeA+BWafWZlVvsVmuI6cF8ab4VLT2NbFptz5qIpYilIpZPJwZP9D6Hgu4QLHjVRqpAruOwQoVagPIwMlXOa6bIolX5jKTnOHjeiVlV8nhPviYhpv89D0rjib4Vuwt4xznwf7zgljEW5ZdtYnF3IOZWfptbXsDPdOpS63stbEynBN5B5I3FkKdSb7Mlxez+XoFPVYK377vflbo5nluT8sQt0UGSqCe2FvJgiX464qeUSUZoFcyvVtpYZaP4MGqlwB0hXTpDt13pQX2Pf8Hd4bbWJOCDcun+LbPNxJk0QoRvpt6eVCYb/F6peT4ZZmZSflJH2SnaXX16CE21CWe05XVsfJcBtqS5t4z8NtK7OW4brsvDCLbnF+C/aetbEclzkdxNoQg43r0/kJV0OfwKSQcvxZGaAZacP/AeFbW8weOAAAIABJREFUuH9hCu25Y19gn7MQ0wv3+XVcuMspw91T/ibyNQ8Pj39dR/fgarWc+zzYxbeG7KwLO3MOZefr9TXwPjMowzXoMnVlRiTBbaqtbbDoFD1LoQ7iXp7h6j4FxGup2dYCM25X9w+7UIflG6oohUKpoQl3ygKc44JEl53lqtZoZnqj74gol693Adwzwm+KH8W5nBPFftMLsOCwcFfFFRN+GxfH89dryzrSby9rMXfqtUK4N/T6/OzsQ9mHDmVn6/WGBt4M90aNTqerrMGtmoH7kbH6FLhnKZjqHGFcp/oUUK+lptHMptx7IkQ/QkMF9u4gGS7DuODRdJmc9RWlKjQAoWTcmIcSLpOG+U0RE547ZSx7B1ovHD3qsHBXk4Ybt45vd+ezlfmAcYn4lsO3vD0KBZdBpnUo+1A2QIJSlHBZGW5dmU6nM5iwfQoNhpo2pA+XczVWZsVve75PgZnt7TDjdps7d6FuKDw04ydMmDBhOvc1UZTjahZoIuVwFjEvYIHYG/eOqHYRvrXcMUI8xY9QMNQnRm5A43Ez15XDwo37EvhtUfFa5Ct0/dhcWlB4uZQ+lZHXcauvV5fo9VnZlw5lZ18CgjXYkuHWVOp0lXXoihnptzUG5qZmKN0yp9c02rbXrAf7FOjZbmQxLjjpGdGPlQEIdwLYM8l6XZRMzjCwgnV2uVKl8RboEoMGyrfQvdGmTCE8d+wQds9Cj4e5jgp3CTDcJ0WCncOrLxWS3WC8hFuiBw6bA1ICPV+PAibD/a+aylxdHX+XAsgSOJrFnRfWVNvnfQpMj4LxW8hvwdglRryQGGtU6HNMhjt9gUZDES81wuURls9UcKKrUmoWBKDvgyddHONS90abSZDu1LHsZMFzBSIH15ajwo1bCS7LcHEC/KIvkTuZkaMkHzhsTlHOpex89ApN2G9vl+oydXUCfbimWi3SiYucFUbtUxdkXIzf2rD3AZmCfQr0mpnxCey3Zv7TQcK84N5ay4AzXIXldOhwWSRMvBLopNEATWigYJcYTAsI3zKcKxYTme6kF9g9C4v4ENJF5ahwg0B3grUlkhfzoHNwiTvr6PX5ly41NzdnZ+dikwXhDPf27a9rKjMrbwn24ZaWmhDVonxLVGOdo4zb2d4O9kFiZouxsw1xWma9F+e3hOd2mWjFUo7bvQi/OuYZq1YrkWcXTCEMF8ppKd9VyiNYryNOuqH8V7lAE8Dj64hyBRLdUDLT5bTm9vQBIQ4Kd2X2yfxY6wBefI88J6wmX38hvSg9PSk/67yuRDhR4KyZNUB+W5ObWYbJcOHRyGQJ1giXOI3Jxgy3zWg0GivBTt6rVzPBTl6whdcAnmT9l3yigtgcebXi6tVzV6/eNbYQ99fh61NgPNfYBfttl9n8AHs6iCpocVBsUBiSBwDCHc9+bnqEQhUaLuPsqFBImV4HeahGhT9lATdwjEueoysmTZfTs9DDsOCgcINO5i+24TDfJbu+/CQ9Pb05J1dvAI6bm6UXWm3A8K2lD7dE91FZA5Lhcgi3rqwNo1kM35J9CugOM3iYjMb6c+eAXCuAMNswO3yFrsTAeAic+W7F1atXKq5erW9vaUe8lu4L6zRSGS7tuOZyHM3OV0cp1epYdRyGcCchvhkul3McV7MgXEK7bFSUZkEA0sfAmy7gGJe5N5pfoHgSO1kYwZs3uaQcE+6Sdfnr/oE8a6lry5Z5xsbGxsYVX9IyiUKV/pAWyRRQ/SJrZsBrb+hO66rwfbiQ997i2C1/hkvRQ1sdhm9NNYar566eu3Ku3mA0ueC+ZgzjEhrOvHruXH0LZ28vGAY243aZsSc9h6kJzg2K9WL1GhCEi7srZbgskns6jkpC34vSnyRdWwfCtxDnhk6d4jduFItyR/domOuYcOOKiv+GPAn2+z63cl3s4pXFbyxbRv51+6mcznCrsrJK7T5HgSJcfaa+xPr9Hu6U1vHsNcMzLvBfA8th6yt05yoq6msx9+5FFMqjU4RrcZlCJ7E7sqKisr7Tsl5mbIcyBWKUY1hWHQ3+G7gu1is2jnDkqECUcKGhlEeolAoF+1mpnPwYqFqgUQYi38NPugjfwvcA9ps0ltMmZuXa3alyTLivcv4yLSsuXhcbWwwEe439ldUXyDs9dGSVIm5rg99W375Zkplb1oD04aKMW6M34faa8Q/gr+DekeBDRS44esZgbMQzrrPnKeD7FLpa2u8C/da3AKftMrAZ19yFo1zaZwPUXkFEvrAYrLFNAoSL7dMlGFc8Xc6+a5o8hNRhQKhG7I3wssDAMS7FuaBn4fe9d59fh4T7jODuZ0uXbSkm6rkttMOiVVRaXWDI1ZUiGS6iVQzjXq/T6XU3sH24nNPCqgw1d/gUi8twKc81VdZX6nT1ZUYTxmPt78VFXdaGTKGzBZwRcre+3sjyW3O3OV2BuiAzROpY9WJl3GJwYhiWcImhkJPuqpqujIISBCV56oJSpVkQiPF1ftJF+BbiXCRZGNmDt4ZwzHEpwRZvQRwWqWXN2iys2/IlCpDfanW5ZdgeBSTDrSmrwu7tFWBcY22lrrIit9bYxmFc/j4Fq46LYwcrvbg04z5uqagAh9s8hjz3nprf+0I1gWFBsUGxciHC1YQqLalCuBK699904uzyeYB2bVw/IwfCt2zOfYFzn9+e2/PrmHAXLVsmdG0G17PybK095ygw3nurTKevEbjfA5Th3qmpaxA4TwHh27aael1lZSW4W6/JiKYJuIEoFKtTxGXxjItmCsQ0mNvb7wLzbaeyhW+LAlAPhIYyVh22QLNgOi/hAseFXBucYU6nw2Ii0VWKF2C7xIQ8FyVcC+dydkMM6TlYcEy4q23F7mfLli9fhzitIOGS2jXodLobDbadpVBXZuI7TQEdJgPQbFldE9WnUNOjfQrwbjOc17L2mhk7iauzlpb6u/V3azu7zOaOmBmI58EjLDZKiHDBqfoRrN0P4nDFdJppI2RiTSA4YcwuyuXnXHDOzZSxI1g9C6N77PrMwRx3nU13ZPtp6ZZrHh6LrfItkuGW6LJu8NzvASHchpoarl45a2YM3zYaKjN11KIDvWZWhyYIjvYpCDMuxmu/Z+8162y38G1X/dWKivbOB0XlahFPR8ECzTwvEfBRfsLlOC5Ju3LadRVS8QJwZx6xsK8jnsub6YaKQydw2sQW9dSt+xwU7pbVyFNoLVlNHJi64hLitfyEC/y2Jje3BHd2DTbDrSu7xXd6DYtxGwyVunNlnBwXtJNDUS5fn0JP9uLC4y57B8939RUVd1vM9/N2qUVYDg0jGBgQ7gT4npGw6shTRtn5llihID12uiQqCqS83ravn9Geix9icehY9jkLPZYsONqrEGv1b9Kq1UuprDe2xJ59ZmWZuluYLgV8hnunstbKeWGEOo2V564auC2MlA/XIV6LDkSfeJ0iLmun5za2d3Nyhe7OlqtX73Y+7ijfpUZ9MVC5QLPAT4BwNQsUkbhkQhxFdumGSpTEu/J6OlvxsOPiR6h4CifMHd1DK7+OCncZuj+SVU+XL2ekvewC4ra8hGs4XlmAPysMR7gmbQOy8wHxW4Mus7IOYV1mrxl5tyj+PgX7MgVk2MG4xOoZmx6I2V5fUdFuftK8K2Y+hkaJe+xMR70W57jQv/3hCuJkMplERdwZAn1fwYHwrYVzJ3LugRbDF5Q6Vw43kq8UzMGWvAH/cotKbCTcmky9UB8ul3ANZQgccAm3VnfO0IjPcKm+MFMd6rB9w7hd3ZWI41Kjk9Bud0feq2EzOF4oRLg4xmVGFLEjQk7cjZKvS4zfc3GMS3DuVM6tfnvo1hAOC3eFwK9n1Rb2/snlKOWi6r1ebTikLxDow+VmuFVUdsvruHdyz+qM+F5cOB2znI3Lx7g93ItrGZ0m1G/p+V1L/dW75u4HHXm7YqIt+oEIF1Ubj+PSSiRYVy6ZTq6huYpzJ3I6c4N65NYQju/yXckbLkOUQNW6Kus9Cobc3FtIl4JAhltjaEB6cWHGrTyXa8DvNePs7cWsmqHuixmITlGXta0XF9r70IJ4LWu01Fe0mLs7H9wvigmjaJUwXESVzIjEMi7DDOFyVZRCOl2zwJu3S4zXc/GUOw9ZPxsW1BMHODou3KX4TuGnL65AFf1GnLVEoSZXZ+A9u4bJcC2duKayKtxeM6pMlRUVhgaEay0DPk+BOv/OecZFyMG2/WbwXjNjJ8ZrIWV3dbeABYrurscdzWp1WCSZ4U5HT8mzOK7Qii5IGJRylUwWatMpC5yBY1wwmXtDULNHVn6dOBCkGLcBYtWyZT8hT3pcu1CCOi3EtyX6XANfHy4uw22oq2lAdprR1VZWqSszWe9TYDy3FuOwvXJmGIdxAcsaEZdFRufd+ooWsLb24F5zs1oyYfx4VJHMiJArkefYQxUwXyOXiecJr59RDo26Lm5O4NzptydOE3NCuG9iLPfN5ajdgirOFyJcXa4WWfdFMgUow72lrcHsNSOrrlKnq2uy0qfAJgeTqYfODLObcc3mWozPQldvwHPN3V3t9fV3O8HjPf7qMWq1iM9vF2giI8KR5ziuq5qhCFdKw207ZYHrudgBbiTFus9vnOuTBacOveMeYfNs+VKeePfFZv4MtySXvibjI1xOhltXU4XZa0bKVldZdkewTwHx29aHdwwo1fZYLy6v35K9uMZOxGGxo6ulvr6lu2vPmjVrvrzXXKQO4yGCyEhrjgt6xMLDlVIF9pQFK56LEi4xkAMcBXeDO1TOCHcJZ1XkxTf4/2KtLOQh3Cpd1i3cOgQmw6UIt7QW6cUl/dZUdk6HScTQ0co5M8wgyLiu6VRAvRbjt2DvGeqzLMaFZv3Vy2vWrFnzd/D4XnmeWoTZ6xPJnAfCP0KngPtCSPGrc1bGPAzjhornjRnBvs/vYpcnC04d7PwG65ezVOgQ9WXpiNMSSi3LNGD6cAXuadZQVoXZa0ZckOUauIq1zrfEaOTmCn3GuN3mlu+wDosd5vysyDV7us3dRLU0l6vDohHGDedXLD0ATWgUUn9+4uD3XBzjakKncFoWhrj8gBCnhPssFnq8BU+3dAXhCPdGVi6ai2H81pLh1mmxZyk01OVe5ew3EzhPgUsPTYYeZ1zEa1HPJdXb3o74LN5vu83dJ9esWROta4eee3Jh16IwmGqtMy7I1KYAPSpnSBzzXBznjhvC7swd7WpYcO4eEKsZy12yhYdumZdeQP22LKsO7cPlEi5zTzPguLU16F4zgLbndCZIscJ7zTCea0C4tq8Y19xtQ65Ad5ABUNjTWV/R2Q3Xg7xdMfPtYVzNAvKMR5W/xDpXYDwXZVzxvJkcyx0WY0Uf9paTd91ZRD1YbrXv8lk5l3ALcnWYdQjuWQpwn8KtshsI24JM9yPdLd41X9xAzsT9ps7mc3EFHRcZ9jOuuaUN9Vp2pkBPLUm4Xd9VXG1nf+3x/fJFYaBzJsoWxl2gERMrvprAiRLk7HJHORfcjQfm3OHcQ+v7VLgeq4l/AZ5h1hyQ+rmD7bdloHcRy7d89zSrqcGdplB3rpJvt5kNGS41DP2AcSkfbW9HvRU7CMNdQxDud3crurn14F55jFoZIY9CdYqOKeQz3oHSaORruIFqF5mcw8SGjXbtNh5n7ywZRHTUCDbcUHWtnEW4uboqHr/ly3DL6tC9Zg0mcPod0qUgeJ4CyritrTrePoXe2G/GPk+BN1fgMO7f14xcs+bv1Oed9RXtCAM/7mjeRa8NC48p5HWZ2FsTIUEc1THO9eXsPxu6CJGEM+WscJcWe7yItCbgayW0elaSWca9HuMhXJpvq0pNyF6zpppKXR2uT4HLt9YYlzhFjGcg+uTRKeqx9vaGUb7a3ojzV3SADHePmfHZlvr6TsR27145fb+5SC3hSXmZoaFOMg8MXKCU2kLFXM9FKddv7Ci25br2PlJO34R68c/LbeyhWHaJIdyyLANypSacKdSUNqBZQmVuLdqnYJVxsfd9MOD51kX7zWxdNyNzBSPitTjG/TsAhVL4uRZyJRie9VeutHd1dtwvUqvnY1Jey5hC9papAjQalVSBeqoDnDvxBXaWO9TLlSu/Tgv3RZv/Hv1UTl+V6XNvcV1WOMO9VVNz5zb3vLA6XZmwYm3pU2ByhTY7GVdApZiBei3quYyn2pYrWAiXqe/u1rdzHfdKC/HgSUdzszpMjiqWGvQdg8G9eEKj5bYkuqh2Off95VruqCBbt4bbUE4L12OdLYBLVPEDgnBLsnRct+X1W7Iv7EZZCbLXzFSZW8daN8OoF8+4GHIAJ45aYVyn9pvZxbjEYUxWM4Una9asGXmSy70t9XdZn9dfyWihH3c+uA86ynj6F+m7nozUaBaEKqRixFPt59xJnDtDDHNlmOu8cHnaGzF1LQ9otCa3BnM9xkO4RIZbUobsNWsqPVeD68Ntwme4vL24vLmC1RQX0SnGZ+1fNyP3mdmSKwDCXdOFQG13y7lOnOPS9bi5PE8tR3WrUVF7IOYR+8+ipLakaNY4l2u5rryPlPPC9Vhn868mtqr6emkmpsNRMMPV1iB7ze7orqJ9YTjt4gaWcVsNbf2FcXlyBba3EoR7EskRwKy4+x3zmGBc7ms6isqDwrhXYKG0E3sTLKyShCOeajfnKjn3P/vDiMUu6yl3gXCf2my5S+9dL82tgvpwUb5FMtzr+v+C/ZZgXFMucpwoTrP2MG5r2zd9u98MPjesHd3tiyVcsxkxXMJ0r3byOi5dD5p3xYj+g8UKNNeSu4k1UpwvW/FcJMvlNOa6EBZcIFyPYpuPKynS4U7QF8pwq/S3kPPCanNN7HQBKTsZl/DgWjzj2uy4CDvYd6YC62QQXK6AZriliJdSs6WihXpcceVKJ/p1cj7u2LUI2vMuplMH6l6pGqn9iS6Xc/3GcCh3aIyrliFcIdwXbbbcdYeweuUj3Nu3a7TIXrM7Oh2+F7fJsT4FfK5gN+MiLus445q7DYjDcgZor9mDIVyqOuvvkg/qr1Sg4S5UT9LLY0QUNTDdOAGU98olVvf+Ip7LodyJXMsd5apbQ7hCuB7P2XKuDUhyV69D/VYgU7hVWoOeppDJvY8k4rgYroXVjGfcVtMTIcbt+TMVYBW3G4UZ9weScHE+Ss27FSTvXrnyHf0c05vDnp0d93fFqAErUImYhtnFo5DYnehyONdvDCfLHeYp1PxqR7lEuKtW2sTcby73KKrGEC5fhnur9gZymkLDFX68dbBPgZMr2NSngNEp6rIOM273d1Z2+xKRwg94wqVNlwDdiitoEwOmHnc0L1LLp9Mu6k1Tg0IqtEcY77lszkUo9zeerkkWXCJcjzesnGtD1IvLPDx+voe4Le/e3prSW9zzwu7U5t5A1s+QEuBbAcYlcwWeXlwbHNeFjIvNFSCf/Ds2w+XMzrMtBONCz3MH63segOU1sq9RxZCvQmp3osvh3DFDOZQ7fDSiDEfKNcJdZQO5XAOXcE+LEK3yEW5pSQNyXpg+14TwrasYt7XRiD6HqpNPp4jLOsO45k7BXOEkYbjdgo7b3f1dRaetjktW+12wNhw2Hz5lQS63IdFFaQGiXM72s2EuOnrUNcL1WG2VXK6tJlpx1mkF/BbKFBq0Jdy9Zg1NNbkmhG/5HJdvIHzbyuw3q+s3jGvu4uYKsDcShLsH9Vju/K6i/eqVCtbzuMF8vbO2u6v7wb2iIrVoNJ2NaSKVUnsTXXbPwpgR7Cx32ChbbjRmtVwk3GfW2oRXrSY3Uq64h+YImAz3VlkVel5YzTkTZq8ZXrF2ZbisXMGxXlxhx8V4rSDjCucKBOG2ma05bnd3Z4VdjtvdTnY6EGvDYf6UciOiZHynOPF5LotzJw7hWu5IgdO7bC4XCddjmXCvzVdM6+NiAcZl/LZKz91r9nXDnZpMtD8M47iCGa4Q47Ya63gZ1wbCdSnjgkOe+TKFH0CGexLxV+y8Isi43Jyh02h5/KB5lzpMQeyi1EREIq4qPIQpd5gX91yDPhSulTPKl/5IP1r5JW+GyzBuTSl8XhjluCU6dK8Zpk/B/l5cZjRizm1E1MmjUWHHRbwW47kcNRlRfdlFuGTZ57jdRvan9/LK1VJVlHKBQirYE4l6rjDljopFBGJ3uUy4PwoFCy9aTvddcZ/Pb+kMt7q0FD0Tt+lWLprfYvzWKuGijGvZa1bb5tS6GTLsYVwuORi+42HcH+h+RrzHsueVKxXfWWNcyHPbjdz3eNxcHqP210yXOHGfiDEjOFnuH0Y6v/LrMuEKnVH+0wpLzvs0rpqPcUm+rS6rwp0Xpqth3dcMW04ybquprt8wLn+uQGe4NjquXZb7XTvyFFheyylfpJZYS3QR9TKcOwmx3GFBNvfC8pXrhLuE//psC/zLLL6H91uKcQu1VbjzFHS5+PMUUP0KMi6GHKC1tIYankzBoTMVnGJcc5cBz7gNIMPdg3orbgLHrb/bgmgfN8jvN2Lfp+Vc54Ny/HnovAOm3LEjuJTrPCy4Trj8B+a+yTrwZlkzH+MCzy0prWbvNSP7FEy6O+ReMwHG5clwbejFFehXQBTKp1PEZZ1k3G7jtzh9kYT7924bHffKlbvf1X+HPM1b3xqxX7nb0t1t1uU171KLEKfl9VyIcrmHLAwbNtrGLgHecqFwr/HdYGUFe0E4TsBvS0sb2H5L9SkcN6FeK+y4+MHTp0BPUx3Wbx3rxnWOcYkoF2Vcak86zhdx8+yV+u72q4j6sYP4fiP+fc6Bjy27Wx7cy4tRR/G5LB/nzhz7x2HcnoUg/nPmbCoXCtejGG+5qzmL0z9fxvAtmeHq67DnhTVVlt22j3AdY9yHbXV9dW4Y5soJmyucZAjXVsft7m7Bt+Riqx3fTNYJSNnc/m69ubuzozmGb887ol46y0Updwg/Wfa6cJ95Ik95eHj8g7vC92IRlnFv37ylLeGcQU4x7q1c5DwFq37rCOO2tRrwjOvAmQrw2bgOMa7Z2IUyLpHh7sF5Im4SjAt6xL5D/wbgRle3+TsD5n26us31j4m1tcxc4vMn98DuNaun3ljuCzF2FJdyh3o6dx8pVwoXf0b5aiRtKMLwbfXtm1XaW9jzwppuV9ZgzgtzaZ8Ckytw7zOJKJRHp6jL2um5qI5wucLfLYRru+N2d9YjX+CtdnyXL2G53eZuXSb1qwI3oihSR3PzXY5yacqdhFLuMC+nkgWXCnfVOrS9cRXaUvHGPdRvr98s0d/CnRfW8HXDrUyEbB0jXKuM+6jBgPFbZzsVHGNcdq5A+l6XPRkuzbjmbnNFJ/p3gGd0onvUiFlPPV/50WPmuccdzUVBYf6I06KcO3PsUG6WO8y5W0O4VLgeb6D2jznF8cU8lHFvlpTizgsjHPcjE3piGFKwPoX5VtBzDX3FuBjPRXMFMlKgP7NendSWM3ss14A8QxaZB5u729NqWb+mjuYYtSiQ13PpLBdnuain2V6uFa5HHFemq1Ape3jEIn57q1SLO7+GzHDfhc9T6JleXEu/wjc9wLg4r7XOuOY2I4dxwaLZyDU4P+SZLQAVzFQqYONo6cS/Vz29j838pweW54lf3+P7YHkN9Vu4awHtWBg6arETyYKLhfsGN577GYUHD48tl7iMW1vCf1+zujLsmbiu7VPA9ysg+uTVKeqyzua4YLcvnnAfP+l40NFx715zc/P95uaiC0XNRUWxRIHHzffBuH+v40HHgw56k+9d7JIYtjp5MohOahubuducm4n+Wp80l6slSpzn8lPuaCfaxFws3Gcr2X+JVmH7dK+Vs/22QFuF81uKcTPvoGwr7Lh8A+FbhHHbHpWhfttXjGs2dloY9/GDjnux6hi1V2xs3BvLli1btkTw2mbVEvCa4vJFQXnl5eU5904/sPGuKGRHJeK33ebuq5bHd3d0sj2XnA+aF8WIkB0TFOeO4zLu0GF/cGIbj4uF67H6r6xP38T/6a5j8W1V2S3WXjNOn8I5GwjX+T4FetS29RvGNbe1d5vvdtwrLy8vio1bsWzF755//vnfIn+YArWsuPhfwVevLVu5fF1sUHn5Fx0PrHpvJ88r7tJWDP4t2M2zK+5xR3mMeg2OcychfblOtYm5WrhP2Qu/WMP18FiRDzFuqfYWel6YhXFrKnkyXH6/dYJxH5nqUMb9775g3AfNzeXlccWrly2laevffvf87/4N+bMUqmXF9KLQctJBXly24rniol15zfcEcgaDGeO3IMuFP39XRz1G/p0wP7kPtl4inIucsQBOE3M4WXC1cD2Ww9vPsJdm4Pk8C9+W1iKnKbD6FMrq6DNxe4Nx2feDQPTJq1PUZR1n3AfNeeXrit9Y9iP7+uC3zz///PMLkT9LoVpWXEz9Y7yM1bu9ZNnq4rjy8uYOxFZBGXmaGyxHNJi7zd/lZiKKt4zvOu7HxITR90Lh68sFyYKjsOBy4XrAwP0jX6t7ucVvSzHn18ApbiVun5mw4/INWxi3tQ5lXJscFxkOMG7n/ebyopVv4G/P+fzzv3v++X9HnhYqi3C/wqwNrVq2pTivvPkel36/I9cYEM9taWF9/uBPvJ5L/I6e3EtfFDODUi7Rl4sw7rChDsOC64X7I7SJh8dwPTzeyCf5trqsBnNeGKtPQddkX5+CM4wLuLaxsS8Y9/GD+/fjvixe/iMuhSFqCTBcrKL5i2Zc8JDnRU+Xrij+Mq753gOoSYEnV+iGWnuBsmv3VaJeSw+aGnKK1JIoMsvFWu5IXo0Il+uFC7U3/jtvo/uqPMJvq7UlvH5LMy5vnwK/4zrMuMBnW8sc2/+ATFsZF+xNjFtZvIJPWlQBwn3+/yJPC5aFcfmuNsh6umzLyrh19+89JlX3hPJgruce4Xz++COdoOcSPRo/dNxXq9XzQccCJssdNtTBNrEeEK7lTqnX+P82lQO+vaWtQs8L4/QpNH3Ek+Ha3Ytra65gaO09xu3suN98r3jLMl6fZWoOMNwfkaeFa0Xxc7QsrF8FvbhoHqTSAAAgAElEQVTiuZX3mzueoHvPUMilPLczNxf1WmowCv7ebH5yr6hZLVJiLdfBW0P0gHA9mJvwvMmHuB4eW/JvV1dpb6DnhXH6FL5uyESdVthv+QamTwHHuI9MJj7GRfQrwA7WGfdxx/2i2GLbLk5++9t/A4Rr7+myK4p/poVr42FzX/28rry5o5bussT3K0Cz4t1Oq55LTbC8NnIoyrleDt3ntyeEy3SUC6xFr7pws1pbjZ4XhjJupr0prpOM++hhZc8z7rcd94vWMfxpreb82/NEcTLcf/lnK9/nsbr4DXoNXhAVOD+b5xYX3X+AaczlNPaSq3t/4j0/HfZccjYvjvUajeS5ArsV+asnhOsRR4eHyFcsFVuirUbOU0D7FJoIx7WjT4GfcfH0gDDuo1ZDDzPug3vNsTxN99j6LaBbMCzC/dd/+d//9Otf/xPu1XBtKWb0auchiW+uLErveIL0K6Dc29W5o8VWz31kXufhsWJd7JiRI2DKdeiAkB4R7t+CyI9ClxvF97AnhnH6FCjH7ZVeXGaYGnuOcR93lAc9Z1828Ox5usjVh3/95//9v35NlFXHhYQrZCI89WNcefMDlsdym8xIZ/0oF3VbHs9tJiTxdNm6RZ7DoTDXAVjoEeF6bCH/mIT+sK5dEPRby4lhOpRuXce4bRjGbWur7SHGfdK8KM6G6zBO/dvztOP+7rf//M//9GtL/QvyWu7PwSJcIRPhryU/l5d3dOLXziwzN7XTRs9tY67EXtwSFDN6FMW5MfbfGqJnhOtBWq7gH1Ys1atwHXNfM/jEMN1/2UO4rF5cBxm31WCFcRGN2sK49/PWrcZ3blippaTf/uY3/+f//zWrrJICJNxrgj8LoXq2fGXefcph7yLYSzpr7bt8K8hczy2C/rV5tvSNxTGjiTBX6DQZfPWQcJcT7Y1CjuuxMh9zBjmGcWtqMXzL77h4vsVnuDyM+8jQ5vzZuCzGfdxcHrfC0Z0qYMXsd7/5/36NlFVSgIS73JnjlJ8tX1nezMe45F6Idytt89zHnHRj1YriGK+RDpwm1kPC9SDuCyT4t3zpfSzhss9TaLrddENnT5+CKxi3td0oTLiIRgUd19icV7zMiXsq/kgaLqJb66QACfe5n5Av2lWkdtFTcWhn7czkSXTNHM+9jwL+0jfiFtm957enhLscmL/woQ+xyBnkN+ndvdCZuA1NZxCvxfqt6xiXmyvY5rjIAKr9oaP5/kr7sZZVP5J8iyj3fyGvRGpLMf0jEPzHz7Z6tmJl8/0nOL8lZm5mpy2eew8bIazqN47rEXdNoFWBqJWXcfc1Q+76cFtX5YpeXHsYlzjRBtGm3Yz7w4Pme8VLnb4l3VMqU/g/HOFaJwWPn4upH8Eqh0J+pFavXgkaIrGe+//aO5feNrIzDfNAnOYtRCOLLGcKyM7b+QP8G9q4mKpTN6K5qgURME6A/gmexaD3ynIgRIbaEzdgps2qntbQzECApJm20yIkNnVpSO1IcSNq0yQzqCpequqcKtaNxSrrPIQ70dWS/Orle77zne/cfetQ0bV67tt+ku6AwMCVlwkXfu5YU7Dc+/BfO7iMi+Baw9VexqQHh4x7cXwcPuNea16LfNNB+NdfGJ6btyk3V1i2gzEXLh9ybswUDgB2IvSvUL/V/nz36Q3Ob22eOwh1uHfOyoQLhNKSqAD6aMK19SlMK7mvXGu4p9H14s4fxyEz7s1lS4hGKxqjn//8F7/4+ceggMTcXN5Vu8r07sSxu4V4ZvrZlMrgBue5P2mTbjAPq+deRvKlrE64UFjiuED4Al9TsN9sdqQifuvsuE4Px/vN7PlW/7N3HiLj/n3QqkcwcxtDDlFuNruRd16kKdMRLTD0VE+dxbkAvt66fot67k+Pd946+u3Mcy8jeR5anXCB0GDdk1VtsKxPYdqJ+0eMTh0065hxPcxUMPXhHp+h2sSr1J5x/3rdr/qvSnpkaE8Lc+3iqxaKERVGEX1B5gZ3UWhd3qCeqz7GzN+xeu7fI5hHvlLhsgLlbrlU/b8xfmu598HoUTjqrqYX15ZxTZ57cRws4972O5Mo7pRxYrjw3GLBouI8JvBSU+FGZLhAsqYf5peda8Rzb353s8xzL6O4L2qFwgUCt6Szo/m1U8Y1JVythrv1yiXjOtRwEZX6ybg/HAXJuFePhIgk4shwY5ZutRfmXQvTxZotNDBKWXtaHkezHMKUN8eTzrTIsHDX73+HPRdhquXeOM2j9cMqhdsW3B0XTJ5i/dZSxdUf+zuhq7hIvnXNuBfHPb8Zt9up86FLX8sp2hpsigWzdjcs2mWMjCtFtUzE/XPC5q69truzs8Rz+0m6dQfHZIJuk5hh+27zFBYZ9/R0+xARKl6xSzLukseiR6F3gmjTNeP2d8vR1CeXUjKEagq1JUvnzUK7w6KibUBQkdzkqOkWrzdKaF1ZPPfum99gHNfkuTchZ+OCVQu3XUGeXKw0XTOu+ayZeuTuuEt6cfEZ1z5TweS5Fwc+Mu71bhVnRiuiOEsKJqzazeYLQwAK2Y8U5WcbxSV5zTuS0xHjj+XKwOK5h789dvJcvZbbQj6Fb1YqXCAvSTPKmXufwmKewpHaw+XblWXc8+MLRJ/4jHs7qPtssA2LVs7NI5+jVLQs1nL5bPYjpfwzLy0N3mi7OPeIF/q3Jn99+1h189zr8Kl7tcIdL7kWCPZxNVxrFXf6OHwZporrO+O+Of7ek+Pe9puy+/e4AgrZLFpC0LOBZbGmO66nngZPyO5ZiGv2TTWGvz1W37p4boAGXBurFS7gsS0Vc4Yd9xqu+bRZt4uoFdFsmF5ca8Y9vzhCVYpk3Ku+wMWwIEMoOLfiFk2hwRBuVJa7NA3B5/1rU0V3x17RNdVyB6Gj1YqFCzruK9o6Lt+a+hTMZ812DmPMuG8ujpfNsTkbVDwenY0e1wabUj5nEa7rO3sGM1oeRVunzebbdP/91tlzK8iH+mTVwoXumzb852gN1+a4i17crYPAvbj+M+4PR28QpzV77s1lNbJVT/QM9SLZP330UTYy4WJGy2MQldbcdW8+PXT03H7YZcGqhQs2Xb9fqoLJuJiEqz9e/f7IwW8jmalwbp1FfnLk4ri3l9VoOr9WhmmlFklU8LyLMWq2pln37u2n106eaz8J4ZuVC5dy/wrrHjLunIudA1fFhjpvZsu4+skzh8dtZ+U7ZKEpLgq7+D4Gn5S9f5bRpDNzXdVW0V0oeDfkT3DlwgW0a8rVp99harh/wd/7YORcDEi+9T1TAakrHDg47ledpuu3lBDma7RIkoK4pCJvRZOukQxuf/sW77lX7hlyKasXbtu1sMD2kYxr61Ow9ijsHIY4b+Yr456f9HB+e9VR1lFH8M+7qXLRim8Qaj6/6bGya0zQ2/71W7zndsL9GFcvXFDF7xROqSzrU7A9uirGcZFMGyTj2u+U7B2gjnvZinGLLBQSB4r5fM6lWdcP0PXfEIuo7F7f/Xj37G+/6WI998pxlKcnYhAu49p/Wf4fh3kKTjc/HG73HDTrmG8DZdw3+p1nlsegxbsuNZPD2H6BsoFcrdRlzQ7piq8DGv6Cwoz2893vv/vu7qetHaznhtv3jUG4QHBbkHIDbMY9xWdc7bGvIh034XtxMRn3vNczO+5ZvymlIyQA0OBwshSrD8qwmqmLmxVuknmAexc8/xL0yFpD2NXqCtf/dotk3B/v+qGaf+IQbnsT++tvMOo49ingMq72p6d2F4qN/rzZ/PH9sWm34UkzjqbFaID4A/F1fTfoQeYhDcYPMhnv9bxG8DE4SrN/c/fT7eND1HPfuteblhCHcIHitgyvYO59QFzW2qdw0t3p4RzX+REk455P6wq62y5rLU4QbR5/YkoyPJbOZBjAZTIZz1sAwYKCDisCvj74+93NVhf13FCbELEId1RxcavJf1pquH9xT7hT9tVvUMVGnXHfTB337LIZ/J8ubt6zvMPz26ZR3nmQeajFN3r6mygrCvbS+wUjbO7whH66eDy5vPrp+23TpJupggd/Cvp54xIumNSdDYtqoX0K/+d878OspnC4vY84bhRzwyx1hZPzN+cXl/WoOrFjoO3Yw8VkdP1RmYypPsk+/IfACf9wvYeBD26MDeNDRaF1+2P38Vu759bd/tolxCPctuDS+1fx0qeAnjd79bL7yqmGi6g0UMZ9c3Hww5tBJUWyHbHSl8grZxgFEdkcbsXNB1qdq5HpIO89B7r7sRujecWLqg7efvfr722eOwjxRBaPcIE8URx/u6pfI724iMuijnt6enq0dWB1XOdHsIx7fjCoJLwlwUKj5mwPM+oZ0z+5kjEWSM2MY8OQywU0S4Em96cqg7e/+dbquTfB7i3RiUm4o+bYccOJHfio4lrpbvdWmXGvNtOzJNMCpZdfskcZU1l9M2N8g+WMU5uhiL3+3usXZHlRevTdp6q1T6zv6GZLiUm4QJqMBadJl31sn8KSjGvQU1UPvbi+5obN/1x1XDerkwbvqe+QMUdccRYbYMZBB+NaiB0Xzna71fiX/d9vWTz3OrjlOnzB0VMHotM+RMVXxrXVwXrP1JVk3OuO47NnEuHwW2UI84hbZXS9Gi9wmQx2T/cTb024eBh0cdDo/+HTG7PnBt89i024nAAYBwcTrlxPmzn77aw01nXMtjNF4zKuw2xc489u8hsXTXiRragPB2lOvZV/YOjVeBvMZFCVaYurMDvcHO6DpUe/vjZ5bvDRjbEJFzQZQOH3TN99HizhzpV5qHaPzpf24i7xXFO+ve0L8Z7bDYXIeWqBeai56tio4gJQ0UxkYhIuJh+L/CfI67zj8EW1ld3fL/bP/hrYcuMTbqMOAIN/6qngMi62TwHnuFpNbE9Vj50zLjY9OGbcs8smzn4Sysect+sdjL0yJZPZ1F5SNkV9TebmuO1QQ3VF7gXyuulXUu3Pe3R/DHxqMj7hAqWtfTu4H0b9a0+9uNiMO/PdY7V7EEXGPR9cys610KTBcF4H2jGZTBmUM7CTkQHTfKDvCkNt9xcYqkYctx1ucpPLXS2i0vpu5rmBh9rEKFxGW0KKNcyPQ5sh5i3hOjiu/jjqqt+i+dZnxr2qlzFfYEJpSJyzPOzwDzKZTQ6IzUwmUzWSEGeqKtj7G6hwQ8wdgsIUtjWY1XL72Ofg5cQoXP1aCHynaCtcxp2fN9vb7p6Eybi3/RAl8ZgRWTlslXk8Ey6PlMNguHkR4rIvjd69NTz3NuARnjiFK+pV7rGENi9VvPbiOjvuNNse7aj7FwEz7lm/in5pyWTE8DgD8MvD6Y6ZkrHNymqgkdcPo+UnfYR/XP5N91yX3WY34hQumBi6gMivo/CV115cRLNor8Lx45e9IBlXraRFtiInR1P1kDNGxNy0bfnWQv4kJA/P/3DzkTYw5G6ApGtPxCrc0TSJM/ap3fCJxyqus+OaH73u9rdoxrVhzbfn+50Qm5txMuIm2AppIAzFytYd3xeTkDVsb08GDN35w493PwYcahOrcMHs/qSxveXmSYCEi8m4iz8H29vHp9799qJVTcWa7D0lRzamGehTCTNVpZKpmn8TGL8neu0wHv1arAud27u7YBMW4hUumEWpEWcNUc0z55lh3vwWmalwsd/d6r72lnEvW2lICSNODnUlLxZ2olhzhxS2iE15blYclwXhs7vbQMd9YxZubZ6lKMteRPmpr15cNONi+xQuDrpGmcHdc7/YTcEBhxGUIHa2aLSIodd8op8eD7nKdW4CLc9iFi5YNL23JfZX81dz/egSruVxsdd9pn574pJxb+tKTEPwA/MJxXEcE8PJ+BEMfQVlG7vF5EijzgmPgvhG3MLlTGU7uPgex/VwVVxrn4K1jnux/436TfcYdVr9MWgmu1X8PcNyHGRiWTiKXOjE5HufWKzzgerRcQvX2IWYwkjz/RXl2ncv7pKMa6mS9Q7UbvfgDMm415eTJC/KxHc8B6m4vkLWZZvWI23/n2IopCHjAsCbG9nG7GyNXBugGTdADddhFaZ78HG32/3m4MRcS3iS4C4wkYIcZONrrmzD8Be9fyL5/XpHDVh+EMTnYxcuUCxaEaemy7ZWk3Ft581Ous+evTw4NjLuVSWpJ3PEBic7TEdYESMYwXSpdtmPbiFXFqq0oATzjviFC23t5KzxWzoIVcV1yrfYOu7ey2fdveOz/toG4bsxYqAkc3E/ETBRbMWNXY5yW96P5cpKlRYmYaok8QsXIOPWeK3RXvDQp+DmuPh8e+qUIHpPOmU5pjWPV0YUV5MlJv7QjWsf8Y+4/LyECOUqTdd5Nry7r0G4FNKCKU7giHsSRZ8CqlB8b9jZ7kSPkZIsw3YCJjCKFC+vRbNa3QLKUfy9DP54yxQRSjRNN/nwOXrKGoQLBLRu2pD51moSLnamwufVxZOa2OBlWWbjKJPi+BUF5do6f3satUhyicOxLN1labpShxFH9nUIt425b/I909mLK+Ned9BdTbHEQkmuOcw5XAVtyEk1iWdjq3ZhoSTKaWqAL9COPwDelQW6Xi7Dd5H8DTbWIVwwQYWj7fpu7S/JuBjmyvWecftNxzWEyDSgxMsyx7Erk9OYhcbfwKw/o4i1RjRfg22+GKs06bpS41yPQYRjLcIdYywXwP7/drf2V55xr7yMsGszFAshB3UaTNjsKTL655O0T8lx1HqiLMJ7iofRyFac3atOQVlRyoJSgysvi6xFuNgOJLHzl9d/VtXuK9Rpo8u4Z33Bz7/VqM0wFAchx0Fe1zHHUgzDUKKr9ESRYXSp6h/D6R+nyX9dMRrLiOKg63fhHe14WgPWFA0ZoguYlbAe4Q5x10I0Na89+ko9eO0/4XrMuNf94H1gQ5EpMUyDhRDyHD+zYw1e/+/sFRxv/H+WYiiGEZN5acToXWSyfSdVNcWWpbgka7Ae4YIaRkC1r/R8e6Sqe5h8G0HGPas3k94HFhMNPgLZjhvcRHdZ3sdZ48hYk3BBBV29s5/NagrdrcNXSL4NmXF/OP+qhQko95D3bNgjZV+ysKxUq8I6+0HXJdwGxnI7C2/t7nRfR5xxB/U0TQNbGSIbZndX260VqrSirP1nuS7hAhq9G7ZpruIebHXPEJk6KXZ5xt1vYX5T7h8i57t/a8oI8mWhWhUmCUlbaxMuhw6C4J9aenH/rG4dLPFbzxn3Mh1nIVfL+3e1IGM+RFai6UqVj29vxgtrEy7mcupG33be7ER9rL6KIOOedTD7OveNMT/x2ygwnDUYIG9ZP+sTLoVOO+ogqfb10c5O19lxEXfFPS4G1ZRcZLpC2InX6Xg6DaVZrSrlEBeXrJr1CRc0kUUCjetTeL2387KHSNaiUOfzZtqew26qRouvAlHy3pHQ0HdrJ5BK+HiUNQq3jTRy88/xfQr76lbXzXGdz5udnj9tIb8f94sxJ3sr2kLll4ZkkbckkTUKFyA3Gootpz6Fiz1V3Tuy59vlGfesn6KpoStAZOHS3QEm7t3aSFincBuC/emogmTcxaOnafc11nEdH/d7z+EFy3ENt2n4LJwogqIoMe/WRsI6hQv+ZF/sV3toxjVxvNNVD05QxTpl3ME93uIdUdC5N3Pc4AyXjWLrdz2sVbhIYQE+xWZcc1Whu71z2LM5Lj7jnrXubRFs3KjJ+BrWCMKyQtNNBf/m9LBW4c4G5s4RWw4Z1+y7F93tZ3snS/PtVeue2u24USvjZMlyE6FOK4kfOOWN9QpXtLc3Yiq52Cquur3dPXLLuKcpGoofJSIrSfY17wvIKU2aTsxubSSsV7iAt2381s/cMq5lVXbU3VL3L07xc8PONnGe84HzvsRPrKO7RFauV+gqF9nZ2uSwZuGObJYLr5wyLo7X+zv/ofYwGfeKvnedYCKUZdMBMpHTGgyaUUwwSCZrFi7grHNtRn2/5816e6oxA9eUcVvK/drj/RJKMmSme2NUWaCryiTq4+BJY93CBbbpILuo1zo77izP7ne3t9WLWb59/ijZg0MjhpVlidJ/UVlF0M7WwsaHFwxQ1i5c1mq51TMPCRczU+H44OX29jdaxr1PDeMNTua0gY7QOA7ufiveh8XahWuaUa4xGWDyrYPj2ntxTw67L9VdZDTZB8mIgjzP8ZKsb33JKWkwiJD1C5ey9NqIT0LNVLit89pJcvaDLuFqk3MVpSkoilBOV4NBhKxfuMD61N5CvNZHN+5T/XO9YPSJBuudbbQaRArKQp2uC1qDwb1uH0qAcEXLVkHza58zFRZ13LMn5rw8ZnmZY6kPpr4wZPhyhaZpoXzP+zQNEiBcUDf/S0yeom7rLeN+gelNEDlJ5mEj5cs1EUrKwwot2K81vM8kQbhWy+0HzLjPq05ONGY5eQIbSZiD65cxa+zW+ppRfy9IgnCBYC6WX6Jui/VbW2LoL9t0EN9JNW0ObkqKnJ9Afbf2nidZZxIhXGDe+BVuEYW6KlbPuBcdb2MT3rdLmnw5Jrn6HTFwQtN0nYtoAugHSjKEq5jOM04GiN8unalwW/VXFBoyDQglWeL8XYO4WkRWatKVijIJf03ePSAZwgWmcgDso/l2Sca9rAczpzFFQchJsgTXWcAXIScL9cqDijK5T1tfIUmIcCcLyx01Eb9d4rj9MvL5/DFmmAavzQmX+NjOsogcx0m1iaRUNjcr92u3NhISItyxaaNW+U9Eog6K1R+3fbQKFvjLYCjGPO/2HUOVGNHtvKFXPmlTTKOhD3rmIa/9lkz0va84hnd/kCREuEBaFAUmAz+O+7y6mkQ41ubfNziOg8ZD0v7Lc5BtNBh9KjnFiBjaFEUZb2eNKebax2gfq3UWsBTFfDyk9LO1wn1sMIiQpAh32JQ46kv9/8K+e761ZNyBEFOFcyRqY/UZqsGysznkhqLn6DPz9YHkUxqagBlxUb9oQGMSsuRrHFIUFD1+jiI6QzOhJEW4gBMAK3Pa4J9x07vjho638TDSJhhU6aaypiRbyCOvwlPIpUW5iREuqGrOxJY55n3TvZK7yLdnyZ/3MZ4O7y6Hubd2SqlULBaxf2zvOMzncgWTAj3rVnvXlCg3OcJtTDd+odT8HLl5B8P56W0nySFxqN8OXo2wwaCYdSBnfb/SRh6UcguxFjZ8iDHvXeRrJTnCBfO0Cj/b+uNhzzHjzuaGXTUTuvs1hny1QlejP6iYw+s2W7S/W2m4kc3OXizZ3+7KcKPg9ubEkCDhzi+nftE8fX2ws60euTruZwkcnDC9HZxb0dnaQjabKyLYPbKQ3dDMea6/DX8eWsxGkGpWT4KEC5TZT+y5rsyjrrrTdZqM29tN1olIpizQtLZbu9JngVwWo6qczVAL2TwAheLiRfRDXMmlIiwkSbji7KbUyRdTp+3tqY/VfUS6p6df7CZmP1+fhCxMVnlv7ZwCkme1QoDtFRtZi/J867DgK1qsiyQJFyjTxZbWrjDPt/t7O9vqK4vjnj5NxFwwdja8G3nLqtCiq11VGzZDHVrfpYh8wDKG2TSk3EQJd2a5YsdaUzg5UFW1O3feV5/7upA3crRJyGvarc1ns/ZEazfUojVO5LO+61u5DeRVySNRwgWTaV2WRrKBNrFG7epDyc+e/GldtUYKymXt3tq17dailmtPuPrazASaLZaST0NWSJZwx9ObUsv72Cru0Z7aVdXqOq7am01Cltdxb60Je8otWnVZKORz2Vy+MF+bFa3P+wVta8L437yjOotpyArJEi6QDMvlLh17Fb7alTnIxqgeVrtquZyU4+B2y7Ua7rBQKGxk84WFcAsWEeYLQ1DcyIN8CYC847ZEKYBLx07ChAtoPb2KHazjnp6eXuv3+jGcLMFVnzwfNzhteLegJGp/rmBJuUjCBcAqbMvTfqFgvCpf0CtrjpabTUHITZpwJSMHtJCMa/B8cTptDCWZo9pe7+/yw2jWYJDAg4qa5S4WX0jCBSXr2ixneqlkGGkhq63XNrLO/TQbWeRViSNpwp3elKr0sI6LjBlva7MKmej0pe3W0nSVh8k9qGhOuQXUGgtZi+rMGswbqSGnW3axMFN0qWhXfy4Fm2eJEy7UGxXhFS7j9rHLsvcM1CZthlTvCErabi1M/CRkc8rdQAzX2DdbYJbx9A05S+wtbhQKBZtx54lwA6D32jAt1HHPWi53m75gYK0mUUEWbW0o0/q8+RfImxLJopaLbJrN/XROFnnWL1nCrbE/UbJ21uSc429iSJ5wef3EbweJt18sv9t0RFEcL3Oel/8NbRKyMOHeIW9JMgvLzWHKVrY1Fypca513uiFcsGxTpKGQmzzhAkHb9aftfvu06rUZoM1ATpbc212g8s90Vamlc978LOUW0YRr81OccPOWqsT0haLFp4njBuKdZrn8tTXjfu5zXvOoTRlHwBqMObRSUJ7eW5vig4ozy0VLCvqGr+XlHCLcDXPEnTfUWEq3AbaJYyeBwgVlTgu5Fsd9EvAcgTaxhtP1y+pXLX8YZ2uNWi4u4YKCbfMAedafW7J+tmfurRa9ozadPJIoXH1GuSXktmpuz/tLaMBys1oXFAmmokPaA0YtF5dw7Wszi3BLQ5MlD3WBz4u+loRBds4Coh2GrSwc93W1FuzzsPpxcCGVt4O7oqVcTA3XHgSsW76FrNYBmZ+q0ijqZheOuxBu0d6BlkQSKVztplT5epZxzzq+l1DjBqcdBxc+kHtrETTLxdRwkWZcLRnMRajbdClnCHe6i4YVboE02QSlBkFjNif3bNOP+rR7a5VqtVn+sO/0yGOOQgDM2gyAjfn75TaKoJQv6uux2Sl0rHDTsP+QUOGOaQAeTbvBKh63soYsT9MVmof34Xq6oT0STCkgz/KLCu0wn89rbWHFXD4/m7qQX2TcRSEhh/2dSBjJFC6QIKjofntFL9ftiOPplRwHTzBFfAzNI3p2b66dvf/QtCuRioibVOGCh0C5PT19dVVH3mJGux2c/jBvBw8Gpj3G9XT6TNYlk7wLKajiJle4/KQxOBOwJeEAAAH+SURBVH39xPFiB0Zpalctk3nzc4q5ErYF3PW05HAagU1Vs5RMBEmqcEFl1D/tT5BXaxtrxtnadDUYrJyStlzD1gPyqJgXGJZb2gg4sGl9JFa4rFDpW7rBmOluLZmEjCWv1QWwGi25Wmghly/kTXsZJZz2E0hihQvqzdmkcQrW9LO1RLIuDDeKTkNCi+4eOrRMeyRD78LCKNoV4bOemHt61bIPSosilx0HQeNIi24TLFxQX8/w7g8RzJlKPD4m6a6ZBAuXEB1IkcyB9LQhEeESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOESUgkRLiGVEOES0gcA4P8BNckzHBrWdmoAAAAASUVORK5CYII=\"\n",
    "height=\"250\"></center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xvlDMc6ocKuz"
   },
   "source": [
    "As the picture above shows, this vector gradient represents the slope of the *plane* that approximates the surface of $f$ near the point $x_0$. The gradient points in the direction of steepest ascent (\"the uphill direction\"), and its magnitude tells us the steepness of the plane (so zero gradient means the plane is flat).\n",
    "\n",
    "**Note**: the gradient points in the direction of *increase* of the function. For gradient descent (discussed next), we will instead be travelling in the direction of *decrease* of the loss. This means that our updates will simply be based on the negative of the gradient."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "bZyrDzor2uCC"
   },
   "source": [
    "## Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "QAZp7A042xH_"
   },
   "source": [
    "Gradient descent is the simplest of training algorithms commonly used in deep learning. However, it gives excellent results in many cases, and also forms the basis for many other powerful optimization methods - such as Momentum, RMSProp, and Adam - which we will look at later in this practical. Mathematically we can describe gradient descent as follows:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{θ}= \\mathbf{θ} −\\eta \\times \\nabla_\\mathbf{θ} J(\\mathbf{θ})\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathbf{θ}$ are the parameters of the model, $\\eta$ (eta) is the learning rate, $J(\\mathbf{θ})$ is the loss (also called $\\mathcal{L}$), and $\\nabla_\\theta J(\\mathbf{θ})$ is the **gradient** of the loss with respect to the parameters (similar to $\\red{\\nabla f(\\mathbf x_0)}$ in the example above). This equation tells us that to update each of the parameters, we scale the gradient for each parameter by the learning rate and subtract it from the corresponding parameter. Or, in pseudo-code:\n",
    "\n",
    "```\n",
    "for each epoch:\n",
    "  grad = calc_grad(loss_func, data, params)  # Calculate gradient of loss wrt parameters\n",
    "  params = params - learning_rate * grad\n",
    "```\n",
    "\n",
    "You might have noticed that we are sweeping a lot of critical details under the rug here! Firstly we are assuming that we can easily calculate the gradients using some `calc_grad` function, and secondly, we are ignoring the issue of batch size (the number of training data points we use to calculate an estimate of the gradient). \n",
    "\n",
    "Luckily for us, TensorFlow addresses the first detail thanks to **automatic differentiation** (AD). With AD, calculating the gradients is about as simple as calling a `calc_grad` function, which means that we do not need to worry about the details of *how* to calculate the gradients. We don't need to think much about implementing derivatives for each of our operations, or the backpropagation algorithm, for example. If you want to know more about how this all works, you should check out the [Build your own TensorFlow](https://colab.research.google.com/drive/14GeXkFd5pQKKNIJ7BMswP0ihlYg5nIbS#forceEdit=true&offline=true&sandboxMode=true) tutorial.\n",
    "\n",
    "In practice, batch size is simply a hyper-parameter that we can tune. However, there are three cases that are worth knowing about:\n",
    "\n",
    "1.   Using a batch size of $n$, where $n$ is the number of training examples is known as **Batch Gradient Descent (BGD)**. In this case, all of the data is used to calculate the gradient at each step. This results in the *most accurate* estimate of the gradient. If the learning rate is not too high, BGD is guaranteed to converge to:\n",
    "    * the global minimum for [convex](https://en.wikipedia.org/wiki/Convex_function) optimisation surfaces\n",
    "    * a local minimum for non-convex surfaces (provided that there are no [saddle points](https://en.wikipedia.org/wiki/Saddle_point)).\n",
    "    \n",
    "    However, some downsides of BGD are that it is **not compatible with online learning**, where we get new examples during training, and that it can be slow for large training datasets.\n",
    "    \n",
    "  **Exercise:** what will happen if we use BGD with too high a learning rate on a convex optimisation surface? Will it diverge or find a local minimum? *Hint:* try drawing a diagram.\n",
    "\n",
    "\n",
    "2.   Using a batch size of 1, called **Stochastic Gradient Descent (SGD)**. In this case, only a single data point is used to calculate an estimate of the gradient. Thus, the estimate is very noisy, and we are not guaranteed to find a minimum (local or global). However, SGD still performs very well in practice and allows for online learning. It also turns out that having noisy estimates of the gradient acts as a form of **regularisation** which can prevent over-fitting. Finally, because we are performing gradient descent with one example at a time, we need much less memory. Not having enough memory can be a significant issue when using BGD.\n",
    "\n",
    "\n",
    "3.   Using a batch size of $m < n$ is called **Mini-batch Gradient Descent** and is a compromise between batch and stochastic gradient descent. We use $m$ examples to calculate an estimate of the gradient. Thus, we still have *some* noise in the gradient estimate, we can tune the batch size to make good use of memory, and the variance of the gradient estimate is greatly reduced — which leads to better convergence to local or global minima.\n",
    "\n",
    "In deep learning, we almost always use mini-batch gradient descent. However, it is often referred to simply as SGD.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Kl8cnkgCSN8_"
   },
   "source": [
    "### Optional extra reading: choosing a batch size\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "n5sBN3xTqKh-"
   },
   "source": [
    "You might be wondering how one chooses which batch size to use for a given problem.\n",
    "\n",
    "One approach is to choose as big a batch size as possible. The reason we might want to do this is that a larger batch size means that our model will train more quickly. This speedup is because modern computer hardware, especially GPUs, are designed with parallelism in mind. In other words, by having a larger batch size, we can take better advantage of our hardware. So in practice, we often choose the largest batch size that will fit into memory.\n",
    "\n",
    "Because the gradient estimation becomes better when more samples are used, using a larger batch size allows you to select a higher learning rate - which means that each iteration would make faster progress. On the other hand, choosing too large a learning rate can be computationally burdensome, and even lead to not converging to a minimum if the descent direction is one of very high curvature. So, in practice we still have to [tune](https://arxiv.org/pdf/1803.09820.pdf) the combination of batch size and learning rate for different tasks and models. \n",
    "\n",
    "[Here](https://arxiv.org/abs/1606.04838) is an excellent reference on optimisation for machine learning that includes a thoughtful comparison of stochastic and full batch gradient methods. For a more advanced read, [this blogpost](https://blog.janestreet.com/does-batch-size-matter/) discusses batch size choice in the context of  recent Bayesian interpretations of SGD.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Y_RP7DKilL0s"
   },
   "source": [
    "Let's look at the one dimensional straight-line aproximation example from above and see how the batch size affects the calculated gradient. In the plot below, we see the function $\\green{f(x)}$, the true tangent line $\\blue{df(x)}$, the true gradient $\\blue{\\tilde{f'}(x_0)}$, the calculated tangent line $\\red{\\tilde{df}(x)}$ and the calculated gradient $\\red{\\tilde{f'}(x_0)}$. In this case, we get the calculated tangent line $\\red{\\tilde{df}(x)}$ and gradient $\\red{\\tilde{f'}(x_0)}$ by using *batch_size* number of datapoints to calculate an estimate.\n",
    "\n",
    "Change the values for $x_0$ and *batch_size* and see what you can learn."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "id": "YRzgfU9ZlL0s"
   },
   "outputs": [],
   "source": [
    "#@title Helper functions (RUN ME) (double click to unhide/hide the code)\n",
    "def f(x):\n",
    "  return -np.cos(x)\n",
    "\n",
    "def tangent_f(x):\n",
    "  return np.sin(x)\n",
    "\n",
    "def df(x, x_0):\n",
    "  return tangent_f(x_0) * (x - x_0) + f(x_0)\n",
    "\n",
    "def perpindicular_unit_f(x_0):\n",
    "  slope_f = tangent_f(x_0)\n",
    "  y_0 = f(x_0)\n",
    "  \n",
    "  x_1 = slope_f / np.sqrt(2) + x_0\n",
    "  y_1 = -x_1 / slope_f + y_0 + x_0 / slope_f\n",
    "  \n",
    "  return [[x_0, x_1], [y_0, y_1]]\n",
    "\n",
    "def noisy_df(x, x_0, noisy_x_0):\n",
    "  return tangent_f(noisy_x_0) * (x - x_0) + f(x_0)\n",
    "\n",
    "def noisy_perpindicular_unit_f(x_0, noisy_x_0):\n",
    "  slope_f = tangent_f(noisy_x_0)\n",
    "  y_0 = f(x_0)\n",
    "  \n",
    "  x_1 = slope_f / np.sqrt(2) + x_0\n",
    "  y_1 = -x_1 / slope_f + y_0 + x_0 / slope_f\n",
    "  \n",
    "  return [[x_0, x_1], [y_0, y_1]]\n",
    "\n",
    "def interactive_noisy_gradient_visual(x_0, noisy_x_0):\n",
    "  # change the fontsize for better visibility\n",
    "  init_size = plt.rcParams[\"font.size\"] # store initial font size\n",
    "  plt.rcParams.update({'font.size': 22}) # update the size\n",
    "  \n",
    "  plt.figure(figsize=(12, 8))\n",
    "  \n",
    "  x = np.linspace(-np.pi, 2 * np.pi)\n",
    "  f_x = f(x)\n",
    "  y_0 = f(x_0)\n",
    "  \n",
    "  # plot f(x)\n",
    "  plt.plot(x, f_x, label=r\"$f(x)$\", color=\"green\")\n",
    "  \n",
    "  # add a point showing where x_0 falls on f(x)\n",
    "  plt.plot(x_0, f(x_0), marker=\"o\", color=\"black\")\n",
    "  \n",
    "  # plot the tangent line to f(x) at x_0\n",
    "  plt.plot(x, df(x, x_0), linestyle=\"--\", color=\"cornflowerblue\", label=r\"$df(x)$\")\n",
    "  \n",
    "  # drop a vertical line from x_0\n",
    "  plt.plot([x_0, x_0], [f(x_0), -3.1], color=\"silver\")\n",
    "  \n",
    "  # plot the noisy tangent line to f(x) at x_0\n",
    "  plt.plot(x, noisy_df(x, x_0, noisy_x_0), linestyle=\"--\", color=\"red\", label=r\"$\\widetilde{df}(x)$\")\n",
    "  \n",
    "  # plot the normal vector to the tangent\n",
    "  [[x_0, x_1], [y_0, y_1]] = perpindicular_unit_f(x_0)\n",
    "  plt.plot([x_0, x_1], [y_0, y_1], color=\"dimgray\")\n",
    "  \n",
    "  # plot the positive direction of change vector\n",
    "  dx = x_1 - x_0\n",
    "  dy = 0 # y_1 - y_1\n",
    "  arrow = plt.arrow(\n",
    "      x_0, y_1, dx, dy,\n",
    "      color=\"blue\", label=r\"$f'(x_0)$\",\n",
    "      lw=3, head_width=np.abs(x_1 - x_0)/10, length_includes_head=True\n",
    "  )\n",
    "  plt.plot([x_0, x_1], [y_1, y_1], color=\"blue\", label=r\"$f'(x_0)$\")\n",
    "  \n",
    "  # plot the noisy normal vector to the tangent\n",
    "  [[noisy_x_0, x_1], [noisy_y_0, y_1]] = noisy_perpindicular_unit_f(x_0, noisy_x_0)\n",
    "  plt.plot([x_0, x_1], [y_0, y_1], color=\"dimgray\")\n",
    "  \n",
    "  # plot the noisy positive direction of change vector\n",
    "  dx = x_1 - x_0\n",
    "  dy = 0 # y_1 - y_1\n",
    "  arrow = plt.arrow(\n",
    "      x_0, y_1, dx, dy,\n",
    "      color=\"red\", label=r\"$f'(x_0)$\",\n",
    "      lw=3, head_width=np.abs(x_1 - x_0)/10, length_includes_head=True\n",
    "  )\n",
    "  plt.plot([x_0, x_1], [y_1, y_1], color=\"red\", label=r\"$\\widetilde{f'}(x_0)$\")\n",
    "  \n",
    "  \n",
    "  plt.legend(loc=\"upper left\")\n",
    "  plt.xlim(-3.1, 6.2)\n",
    "  plt.ylim(-3.1, 3.1)\n",
    "  plt.xlabel(r\"$x_0$\")\n",
    "  plt.show()\n",
    "  \n",
    "  # reset to initial font size\n",
    "  plt.rcParams.update({'font.size': init_size})\n",
    "\n",
    "def interactive_batch_size_visual(x_0, batch_size):\n",
    "  np.random.seed(0)\n",
    "  noisy_x_0 = x_0 + np.mean(np.random.normal(loc=0, scale=0.5, size=batch_size))\n",
    "  interactive_noisy_gradient_visual(x_0, noisy_x_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "c-tKrCuHlL0z"
   },
   "outputs": [],
   "source": [
    "#@title Double click to unhide/hide the code {run: \"auto\"}\n",
    "x_0 = 2.7 #@param {type:\"slider\", min:-3.1, max:6.2, step:0.01}\n",
    "batch_size = 65 #@param {type:\"slider\", min:1, max:256, step:1}\n",
    "\n",
    "interactive_batch_size_visual(x_0, batch_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "0TOg0zsClL02"
   },
   "source": [
    "**Exercise:** What are the general relationships you notice between the accuracy of the calculated gradient $\\red{\\tilde{f'}(x_0)}$ and (a) the batch size and (b) the funciton $\\green{f(x)}$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kqzkyzvvSUxY"
   },
   "source": [
    "### Implementing SGD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "pF-1lI-S0x7r"
   },
   "source": [
    "SGD is incredibly simple to implement, and as you'll see later can perform very well!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SXEq8950yU4Z"
   },
   "outputs": [],
   "source": [
    "def SGD_update(params, grads, states, hyper_params):\n",
    "  # hyper-param typical values: learning_rate=0.01\n",
    "  # SGD doesn't have any state, however, the algorithms\n",
    "  # we will look at later do!\n",
    "  for param, grad in zip(params, grads):\n",
    "    param.assign_sub(hyper_params['lr'] * grad) # TF function that updates by substracting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DHavR63e1bQ1"
   },
   "source": [
    "And now we can visualize each step of the SGD optimization of Rosenbrock's banana function:\n",
    "\n",
    "**Exercise:** Take a look at the hidden code below and make sure that you understand it. Do you see how the code relates to the pseudo-code and equation above? For example, where do you find `calc_grad`? Once you understand the code, you should hide it again to make it easier to change the sliders and see the results.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "id": "kqnxFxuyEweK"
   },
   "outputs": [],
   "source": [
    "#@title Helper functions (RUN ME) (double click to unhide/hide the code)\n",
    "\n",
    "\n",
    "def optimize_banana(update_func, params, states, hyper_params, add_gradient_noise=False):\n",
    "  # plot the loss surface, minimum value and starting point\n",
    "  X, Y, Z = gen_2d_loss_surface(rosenbrock_banana)\n",
    "  fig, ax = make_contour_plot(X, Y, Z)\n",
    "  ax.plot(1, 1, 'r*', ms=30, label='minimum') \n",
    "  ax.plot(start_x, start_y, 'b*', ms=20, label='start')\n",
    "\n",
    "  for epoch in range(epochs):\n",
    "    with tf.GradientTape() as tape:\n",
    "      # we are trying to minimize the output of the 🍌 func\n",
    "      if add_gradient_noise:\n",
    "        # add random noise to the function parameters, this will result in a noisy gradient        \n",
    "        # we only perturb the linear parameter b so in expectation we get the true gradient\n",
    "        loss = rosenbrock_banana(x, y, b=(20.0+np.random.normal())) \n",
    "      else:\n",
    "        loss = rosenbrock_banana(x, y) \n",
    "      \n",
    "    # calculate the gradients of the loss with respect to the params\n",
    "    grads = tape.gradient(loss, params)\n",
    "    \n",
    "    # save the old x and y values for the plot\n",
    "    old_x = params[0].numpy()\n",
    "    old_y = params[1].numpy()\n",
    "\n",
    "    # update the parameters using SGD\n",
    "    update_func(params, grads, states, hyper_params)\n",
    "\n",
    "    # plot the change in x and y for each update step\n",
    "    ax.annotate('', xy=(x.numpy(), y.numpy()),\n",
    "                xytext=(old_x, old_y),\n",
    "              arrowprops={'arrowstyle': '->', 'color': 'k', 'lw': 1},\n",
    "                   va='center', ha='center')  \n",
    "\n",
    "  ax.plot(x.numpy(), y.numpy(), 'g*', ms=20, label='end')\n",
    "  ax.legend()\n",
    "\n",
    "  fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "id": "veRWE23R4TSi"
   },
   "outputs": [],
   "source": [
    "#@title Double click to unhide/hide the code {run: \"auto\"}\n",
    "start_x = -1 #@param {type:\"slider\", min:-2, max:2, step:0.1}\n",
    "start_y = 0.73334 #@param {type:\"slider\", min:-0.26666, max:1.0666, step:0.1}\n",
    "learning_rate = 0.015 #@param {type:\"slider\", min:0, max:0.02, step:0.0005}\n",
    "epochs = 100 #@param {type:\"slider\", min:1, max:150, step:1}\n",
    "\n",
    "x = tf.Variable(start_x, dtype='float32')  \n",
    "y = tf.Variable(start_y, dtype='float32')\n",
    "params = [x, y]\n",
    "states = []\n",
    "hyper_params = {\"lr\": learning_rate}\n",
    "\n",
    "optimize_banana(SGD_update, params, states, hyper_params, add_gradient_noise=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N3Cz4chcpKum"
   },
   "source": [
    "**Remark:** You might have noticed that we are computing the true gradient of the 🍌 function, which has no stochasticity! To make it look more like an actual case of SGD, we added the option to add some noise to the value of $b$ to obtain a stochastic gradient whose expectation is the actual gradient. This is similar to what happens in machine learning when we take random training samples: we can think of the random data points as the random parameters of this loss function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "8j77_Fe02wqe"
   },
   "source": [
    "**Partner Exercise:** Tweak the starting position (x and y), as well as the learning rate and the number of epochs. See if you can get to the global minimum. What do you notice about the behaviour of SGD in the 🍌 valley?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FMJUFyny4EnB"
   },
   "source": [
    "## SGD with Momentum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AzDIrG6w4yQ9"
   },
   "source": [
    "Hopefully, you've noticed a bit of an issue with SGD. When the gradient oscillates sideways or is very small (i.e. in the 🍌) progress towards the minimum is very slow. One solution to this problem is to add a momentum term to our optimization step:\n",
    "\n",
    "\\begin{align}\n",
    "\\Delta \\mathbf{θ} &= \\gamma \\Delta \\mathbf{θ} + \\eta \\nabla_\\mathbf{θ} J(\\mathbf{θ}) \\\\\n",
    "\\mathbf{θ} &= \\mathbf{θ} − \\Delta \\mathbf{θ}\n",
    "\\end{align}\n",
    "\n",
    "where $\\Delta \\mathbf{θ}$ is the change in parameters $\\mathbf{θ}$ at each step and is made up of a mixture between the gradients at a given step and the change from the previous step. $\\gamma$ (gamma) is called the *momentum* term, and $\\eta$ is called the learning rate, as before.\n",
    "\n",
    "The reason that this method is called *momentum* is that we can compare it to SGD as follows:\n",
    "\n",
    "> Gradient descent is a person walking down a hill. They follow the steepest path downwards; their progress is slow, but steady. Momentum is a heavy ball rolling down the same hill. The added inertia acts both as a smoother and an accelerator, dampening oscillations and causing it to barrel through narrow valleys, small humps and local minima. \n",
    "\n",
    "In other words, **the momentum term speeds up optimisation if the direction of change stays more or less the same and reduces oscillations when the direction of change goes back and forth**.\n",
    "\n",
    "We can also describe momentum with some simple pseudo-code:\n",
    "\n",
    "```\n",
    "change = 0\n",
    "for each epoch:\n",
    "  grad = calc_grad(loss_func, data, params)\n",
    "  change = momentum * change + learning_rate * grad\n",
    "  params = params - change\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "WBGPpJHn4vRW"
   },
   "source": [
    "### Implementing Momentum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "KSJmsRTKg86_"
   },
   "outputs": [],
   "source": [
    "def momentum_update(params, grads, states, hyper_params):\n",
    "  # hyper-param typical values: learning_rate=0.01, momentum=0.9\n",
    "  changes = states['changes']\n",
    "  for param, grad in zip(params, grads):\n",
    "    changes[param].assign(hyper_params['momentum'] * changes[param] +\n",
    "                          hyper_params['lr'] * grad)\n",
    "    param.assign_sub(changes[param])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "id": "iDyzTYPlB31-"
   },
   "outputs": [],
   "source": [
    "#@title Double click to unhide/hide the code{run: \"auto\"}\n",
    "start_x = -1 #@param {type:\"slider\", min:-2, max:2, step:0.1}\n",
    "start_y = 0.73334 #@param {type:\"slider\", min:-0.26666, max:1.0666, step:0.1}\n",
    "learning_rate = 0.01 #@param {type:\"slider\", min:0, max:0.02,step:0.0005}\n",
    "momentum = 0.8 #@param {type:\"slider\", min:0, max:0.99, step:0.01}\n",
    "epochs = 100 #@param {type:\"slider\", min:1, max:150,step:1}\n",
    "\n",
    "x = tf.Variable(start_x, dtype='float32')  \n",
    "y = tf.Variable(start_y, dtype='float32')\n",
    "params = [x, y]\n",
    "changes = {param: tf.Variable(0., dtype='float32') for param in params}\n",
    "hyper_params = {\"lr\": learning_rate, \"momentum\": momentum}\n",
    "states = {\"changes\": changes}\n",
    "\n",
    "optimize_banana(momentum_update, params, states, hyper_params, add_gradient_noise=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N8popkDD1IVo"
   },
   "source": [
    "**Partner Exercise:** Play around with the various parameters and see if you can get to the minimum. Compare the performance of momentum in the 🍌 to that of SGD.\n",
    "\n",
    "**Partner Exercise:** Why do we see oscillations for high enough momentum values?\n",
    "\n",
    "**Partner Exercise:** Do you see any relationship between the amount of momentum and the learning rate?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IujnQlvqFpgN"
   },
   "source": [
    "## RMSProp (Root Mean Square Propagation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "J5McgKZaFpBy"
   },
   "source": [
    "So far, we have chosen a learning rate $\\eta$, and it has been multiplied by the whole gradient. Thus, each element in the gradient vector has had the same learning rate applied to it at every step. However, we can imagine that:\n",
    "\n",
    "1.   Each weight might not need to vary by the same amount.\n",
    "2.   The amount we want to change each parameter might change throughout the optimisation process. \n",
    "\n",
    "**Exercise:** can you think of examples of when these two cases might apply?\n",
    "\n",
    "RMSProp is a method that addresses these issues. It can be described with the following formulae:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{v} &= \\gamma \\mathbf{v} + (1 - \\gamma) (\\nabla_\\mathbf{θ} J(\\mathbf{θ}))^2 \\\\\n",
    "\\mathbf{θ} &= \\mathbf{θ} − \\frac{\\eta}{\\sqrt{\\mathbf{v}} + \\mathbf{\\epsilon}} \\nabla_\\mathbf{θ} J(\\mathbf{θ})\n",
    "\\end{align}\n",
    "\n",
    "where each element of $\\mathbf{v}$ is an estimate of the square of the gradient for a specific parameter, calculated using a rolling average, $\\gamma$ is a forgetting factor for the rolling average, and $\\mathbf{\\epsilon}$ (epsilon) is a small number added for numerical stability. \n",
    "\n",
    "In words, RMSProp is scaling down the learning rate for each gradient by rolling average of the most recent gradients for that parameter. Importantly **each parameter has its own learning rate, which changes over time**.\n",
    "\n",
    "We can also describe RMSPRop using pseudo-code:\n",
    "\n",
    "```\n",
    "average = 0\n",
    "for each epoch:\n",
    "  grad = calc_grad(loss_func, data, params)\n",
    "  average = gamma * average + (1 - gamma) * pow(grad, 2)\n",
    "  params = params - learning_rate / sqrt(average) * grad\n",
    " ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "woMqJYFbdGWO"
   },
   "source": [
    "### Implementing RMSProp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DuPoytCSdztD"
   },
   "outputs": [],
   "source": [
    "def RMSProp_update(params, grads, states, hyper_params):\n",
    "  # hyper-param typical values: learning_rate=0.001, gamma=0.9, eps=1e-8\n",
    "  averages = states['averages']\n",
    "  for param, grad in zip(params, grads):\n",
    "    averages[param].assign(hyper_params['gamma']*averages[param] +\n",
    "                           (1 - hyper_params['gamma'])*tf.math.pow(grad, 2))\n",
    "    param.assign_sub(hyper_params['lr']/(tf.sqrt(averages[param]) + hyper_params['eps'])*grad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "id": "xuqG4shffumO"
   },
   "outputs": [],
   "source": [
    "#@title Double click to unhide/hide the code {run: \"auto\"}\n",
    "start_x = -1 #@param {type:\"slider\", min:-2, max:2, step:0.1}\n",
    "start_y = 0.73334 #@param {type:\"slider\", min:-0.26666, max:1.0666, step:0.1}\n",
    "learning_rate = 0.01 #@param {type:\"slider\", min:0, max:0.02,step:0.0005}\n",
    "gamma = 0.63 #@param {type:\"slider\", min:0.01, max:0.99, step:0.01}\n",
    "epochs = 100 #@param {type:\"slider\", min:1, max:150,step:1}\n",
    "\n",
    "x = tf.Variable(start_x, dtype='float32')  \n",
    "y = tf.Variable(start_y, dtype='float32')\n",
    "params = [x, y]\n",
    "averages = {param: tf.Variable(0., dtype='float32') for param in params}\n",
    "states = {\"averages\": averages}\n",
    "hyper_params = {\"lr\": learning_rate, \"gamma\": gamma, \"eps\": 1e-8}\n",
    "  \n",
    "optimize_banana(RMSProp_update, params, states, hyper_params, add_gradient_noise=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "TwjIfkp16Etl"
   },
   "source": [
    "**Partner Exercise:** Play around with the various parameters and see if you can get to the minimum. Compare the performance of RMSProp with momentum and SGD, particularly in the 🍌."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NRymHjfYDagF"
   },
   "source": [
    "## Adam (Adaptive moment estimation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GOtreVpxFE4X"
   },
   "source": [
    "Adam combines the ideas of momentum and adaptive learning rates that we have explored above. More specifically, in addition to storing the rolling averages of the *squared* gradients and using them to control the learning rate for each parameter, like RMSProp, it also stores the rolling averages of the gradients themselves and uses them like momentum. Mathematically, we can describe Adam as follows:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf m &= β_1 \\mathbf{m} + (1 − β_1)  \\nabla_\\mathbf{θ} J(\\mathbf{θ}) \\\\\n",
    "\\mathbf v &= β_2 \\mathbf{v} + (1 − β_2) (\\nabla_\\mathbf{θ} J(\\mathbf{θ}))^2 \\\\\n",
    "\\mathbf{\\hat{m}} &= \\frac{\\mathbf{m}}{1 − β_1^t} \\\\\n",
    "\\mathbf{\\hat{v}} &= \\frac{\\mathbf{v}}{1 − β_2^t} \\\\\n",
    "\\mathbf{θ} &= \\mathbf{θ} − \\frac{\\eta}{\\sqrt{\\mathbf{\\hat{v}}} + \\mathbf{\\epsilon}} \\mathbf{\\hat{m}}\n",
    "\\end{align}\n",
    "\n",
    "where $\\mathbf{m}$ is a rolling estimate of the gradient, $\\mathbf{v}$ is a rolling estimate of the squared gradient,  $\\mathbf{\\hat{m}}$ and $\\mathbf{\\hat{v}}$ are bias-corrected estimates, and the $\\beta_i$ (beta) are mixing factors. As before, $\\eta$ and $\\epsilon$ are the learning rate and a numerical stability term.\n",
    "\n",
    "The name Adam comes from the fact that $\\mathbf{m}$ and $\\mathbf{v}$ are estimates of the first moment (the mean) and the second moment (the uncentered variance) of the gradient. The reason that we need bias corrected versions of the estimates is that they are initialized to be zero vectors, which means that when the training starts, the estimates are biased towards zero. This problem is especially relevant when $β_1$ and $β_2$ are close to 1.\n",
    "\n",
    "The pseudo-code for Adam is slightly longer than the other methods we've looked at but should be reasonably simple to understand. If something doesn't make sense, go back and look at momentum and RMSProp.\n",
    "\n",
    "```\n",
    "first_moment = 0\n",
    "second_moment = 0\n",
    "t = 0\n",
    "for each epoch:\n",
    "  t = t + 1\n",
    "  grad = calc_grad(loss_func, data, params)\n",
    "  first_moment  = beta1 * first_moment  + (1 - beta1) * grad\n",
    "  second_moment = beta2 * second_moment + (1 - beta2) * pow(grad, 2)\n",
    "  first_moment_unbiased  = first_moment  / (1 - pow(beta1, t))\n",
    "  second_moment_unbiased = second_moment / (1 - pow(beta2, t))\n",
    "  params = params - (learning_rate / sqrt(second_moment_unbiased)) * first_moment_unbiased\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "bgvb-C-0q4If"
   },
   "source": [
    "### Implementing Adam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IsCRgDIgeMFI"
   },
   "source": [
    "Implementing Adam is left as an exercise, we've set up the scafolding for you but you'll have to get your hands dirty! So skip this section for now, and come back later!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8Mnm5L8ODbD7"
   },
   "outputs": [],
   "source": [
    "def Adam_update(params, grads, states, hyper_params):\n",
    "  # hyper-param typical values: learning_rate=0.001, beta1=0.9, beta2=0.999, eps=1e-8\n",
    "  t = states[\"t\"]\n",
    "  t.assign_add(1.0)  \n",
    "  first_moments = states[\"first_moments\"]\n",
    "  second_moments = states[\"second_moments\"]\n",
    "  \n",
    "  for param, grad in zip(params, grads):\n",
    "    # come back here later!\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "id": "LCPLzIHDf81t"
   },
   "outputs": [],
   "source": [
    "#@title Double click to unhide/hide the code {run: \"auto\"}\n",
    "start_x = -1 #@param {type:\"slider\", min:-2, max:2, step:0.1}\n",
    "start_y = 0.73334 #@param {type:\"slider\", min:-0.26666, max:1.0666, step:0.1}\n",
    "learning_rate = 0.01 #@param {type:\"slider\", min:0, max:0.02,step:0.0005}\n",
    "beta1 = 0.8 #@param {type:\"slider\", min:0.01, max:0.99, step:0.01}\n",
    "beta2 = 0.625 #@param {type:\"slider\", min:0.01, max:0.999, step:0.001}\n",
    "epochs = 100 #@param {type:\"slider\", min:1, max:150,step:1}\n",
    "\n",
    "x = tf.Variable(start_x, dtype='float32')  \n",
    "y = tf.Variable(start_y, dtype='float32')\n",
    "params = [x, y]\n",
    "first_moments = {param: tf.Variable(0., dtype='float32')\n",
    "                 for param in params}\n",
    "second_moments = {param: tf.Variable(0., dtype='float32')\n",
    "                                     for param in params}\n",
    "t = tf.Variable(0., dtype='float32')\n",
    "states = {\"first_moments\": first_moments, \n",
    "          \"second_moments\": second_moments, \"t\": t}\n",
    "hyper_params = {\"lr\": learning_rate, \"beta1\": beta1, \"beta2\": beta2, \"eps\": 1e-8}\n",
    "\n",
    "optimize_banana(Adam_update, params, states, hyper_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GbDpI0ehaZC2"
   },
   "source": [
    "## Learning Rate Decay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Deoi1TZcdF9-"
   },
   "source": [
    "One of the advantages of methods such as Adam and RMSProp is that the effective learning rate for each parameter is no longer a constant during training. This behaviour is desirable because we often want to reduce the learning rate as we near the global minimum so that we do not overshoot it. Another simple strategy for solving this problem is *learning rate decay*. With learning rate decay, we progressively reduce the learning rate during the training process. For example, we might decay our learning rate as follows:\n",
    "\n",
    "\\begin{equation}\n",
    "  \\eta = \\eta \\times \\frac{1}{1 + \\delta \\times t}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\eta$ is the learning rate, $\\delta$ (delta) is the decay rate, and $t$ is the number of the current training epoch. This method will give us a learning rate that looks something like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "hmzu-nI5iR1j"
   },
   "outputs": [],
   "source": [
    "initial_learning_rate = 0.01\n",
    "epochs = 100\n",
    "decay_rate = initial_learning_rate/epochs\n",
    "learning_rates = [initial_learning_rate]\n",
    "\n",
    "for t in range(epochs):\n",
    "  previous_learning_rate = learning_rates[t]\n",
    "  learning_rate_at_t = previous_learning_rate * 1/(1 + decay_rate * t)\n",
    "  learning_rates.append(learning_rate_at_t)\n",
    "  \n",
    "plt.plot(learning_rates)\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.ylabel(\"LR\")\n",
    "plt.title(\"LR vs Epoch\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aVZ-scjAj7Zr"
   },
   "source": [
    "In practice, while we do often use the simple decay scheme shown above, it is also common to use more complicated *learning rate schedules* such as exponential decay and step decay:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FRY6xZn_mVHG"
   },
   "outputs": [],
   "source": [
    "# exponential decay\n",
    "initial_learning_rate = 0.01\n",
    "epochs = 100\n",
    "decay_rate = 0.01\n",
    "learning_rates = []\n",
    "\n",
    "for t in range(epochs):\n",
    "  learning_rate_at_t = initial_learning_rate * np.exp(-decay_rate * t)\n",
    "  learning_rates.append(learning_rate_at_t)\n",
    "  \n",
    "plt.plot(learning_rates)\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.ylabel(\"LR\")\n",
    "plt.title(\"LR vs Epoch\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "U4_-JDBknfGC"
   },
   "outputs": [],
   "source": [
    "# step decay\n",
    "initial_learning_rate = 0.01\n",
    "epochs = 100\n",
    "epochs_wait = 10\n",
    "decay_rate = 0.5\n",
    "learning_rates = []\n",
    "\n",
    "for t in range(epochs):\n",
    "  exponent = np.floor((1 + t) / epochs_wait) # this part gives us the \"steps\"\n",
    "  learning_rate_at_t = initial_learning_rate * np.power(decay_rate, exponent)\n",
    "  learning_rates.append(learning_rate_at_t)\n",
    "  \n",
    "plt.plot(learning_rates)\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.ylabel(\"LR\")\n",
    "plt.title(\"LR vs Epoch\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UKTj5p9UprUm"
   },
   "source": [
    "Any of these learning rate decay methods are compatible with the optimisation methods described above. The choices of whether or not to use learning rate decay and if so, which method to use, are hyper-parameters that we can tune.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HRtm9JmFD3Ks"
   },
   "source": [
    "## Putting it all into practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "klDdVCbiP0Cn"
   },
   "source": [
    "Putting all of this into practice is very simple! Keras provides us with a high-level API that makes using any of the optimization methods or learning rate schedules as easy as changing a single line of code. Of course, if you are defining a custom training loop using `tf.GradientTape` then the code we've used above can easily be converted to work for any model and dataset. \n",
    "\n",
    "One of the advantages of using Keras is that it provides us with reasonable default values for all of the hyper-parameters of the optimization algorithms. Let's use Keras to train a simple MLP on FashionMNIST so that we can compare the optimization algorithms we've looked at in a more realistic setting.\n",
    "\n",
    "As a quick reminder, FashionMNIST contains 28x28 grayscale images from 10 different types of clothing. Let's take a quick look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "UQxsa6H0SCvo"
   },
   "outputs": [],
   "source": [
    "fashion_mnist = tf.keras.datasets.fashion_mnist\n",
    "\n",
    "(train_images, train_labels), (test_images, test_labels) = fashion_mnist.load_data()\n",
    "\n",
    "text_labels = ['T-shirt/top', 'Trouser', 'Pullover', 'Dress', 'Coat', 'Sandal', 'Shirt', 'Sneaker', 'Bag', 'Ankle boot']\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "for i in range(25):\n",
    "    plt.subplot(5,5,i+1)\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "    plt.grid(False)\n",
    "\n",
    "    img_index = np.random.randint(0, 50000)\n",
    "    plt.imshow(train_images[img_index], cmap=\"gray_r\")\n",
    "    plt.xlabel(text_labels[train_labels[img_index]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_SNA2x5ETHhF"
   },
   "source": [
    "Before we train on the data, we want to do a little pre-processing. We won't go into detail about how and why we are doing the pre-processing, but if you want to know more, you should take a look at the *Deep Feedforward Networks* practical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6YPa_8xcTGn0"
   },
   "outputs": [],
   "source": [
    "batch_size = 128\n",
    "\n",
    "train_ds = tf.data.Dataset.from_tensor_slices((train_images, train_labels))\n",
    "# Divide image values and cast to float so that they end up as a floating point number between 0 and 1\n",
    "train_ds = train_ds.map(lambda x, y: (tf.cast(x, tf.float32) / 255.0, tf.cast(y, tf.int32)))\n",
    "# Shuffle the examples.\n",
    "train_ds = train_ds.shuffle(buffer_size=batch_size * 10)\n",
    "# Now \"chunk\" the examples into batches\n",
    "train_ds = train_ds.batch(batch_size)\n",
    "\n",
    "test_ds = tf.data.Dataset.from_tensor_slices((test_images, test_labels))\n",
    "test_ds = test_ds.map(lambda x, y: (tf.cast(x, tf.float32) / 255.0, tf.cast(y, tf.int32)))\n",
    "test_ds = test_ds.batch(batch_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zLKzrbhRTv2V"
   },
   "source": [
    "Now let's define a simple MLP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "AI1Slx66S78W"
   },
   "outputs": [],
   "source": [
    "def build_mlp():\n",
    "  model = tf.keras.Sequential([\n",
    "    tf.keras.layers.Flatten(input_shape=(28, 28)),\n",
    "    tf.keras.layers.Dense(128, activation='relu'),\n",
    "    tf.keras.layers.Dense(10, activation='softmax')\n",
    "  ])\n",
    "  \n",
    "  return model\n",
    "\n",
    "model = build_mlp()\n",
    "\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UuKDP3YYUXDW"
   },
   "source": [
    "And finally, lets train the MLP using a few different optimisers and compare the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "id": "YI1OL9Orp1Hl"
   },
   "outputs": [],
   "source": [
    "#@title Helper functions (double click to unhide/hide the code)\n",
    "\n",
    "def make_loss_plots(losses):\n",
    "  plt.close()\n",
    "  for label, loss_vals in losses.items():\n",
    "    plt.plot(loss_vals, label=label)\n",
    "  plt.xlabel(\"Epoch\")\n",
    "  plt.ylabel(\"Loss\")\n",
    "  plt.title(\"Loss vs Epoch\")\n",
    "  plt.legend()\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "_-4f4X8XVH2O"
   },
   "outputs": [],
   "source": [
    "losses = {}\n",
    "tf.random.set_seed(0)\n",
    "\n",
    "# SGD\n",
    "model = build_mlp()\n",
    "model.compile(optimizer='sgd', \n",
    "              # we can use the string shortcut if we \n",
    "              # don't want to change the hyper-parameters\n",
    "              loss='sparse_categorical_crossentropy',\n",
    "              metrics=['accuracy'])\n",
    "sgd_hist = model.fit(train_ds, epochs=5,\n",
    "          validation_data=test_ds)\n",
    "\n",
    "losses['SGD'] = sgd_hist.history[\"val_loss\"]\n",
    "make_loss_plots(losses)\n",
    "\n",
    "\n",
    "\n",
    "# Momentum\n",
    "model = build_mlp()\n",
    "model.compile(optimizer=tf.keras.optimizers.SGD(lr=0.01, momentum=0.9),\n",
    "              loss='sparse_categorical_crossentropy',\n",
    "              metrics=['accuracy'])\n",
    "momentum_hist = model.fit(train_ds, epochs=5,\n",
    "          validation_data=test_ds)\n",
    "\n",
    "losses['Momentum'] = momentum_hist.history[\"val_loss\"]\n",
    "make_loss_plots(losses)\n",
    "\n",
    "# RMSProp\n",
    "model = build_mlp()\n",
    "model.compile(optimizer=tf.keras.optimizers.RMSprop(lr=0.001, rho=0.9),\n",
    "              # rho is the symbol used for the forget factor in Keras\n",
    "              # we used gamma (γ) in this practical\n",
    "              loss='sparse_categorical_crossentropy',\n",
    "              metrics=['accuracy'])\n",
    "rmsprop_hist = model.fit(train_ds, epochs=5,\n",
    "          validation_data=test_ds)\n",
    "\n",
    "losses['RMSProp'] = rmsprop_hist.history[\"val_loss\"]\n",
    "make_loss_plots(losses)\n",
    "\n",
    "# Adam\n",
    "# UNCOMMENT THIS ONCE YOU HAVE IMPLEMENTED ADAM!\n",
    "# model = build_mlp()\n",
    "# model.compile(optimizer=tf.keras.optimizers.Adam(lr=0.001, beta_1=0.9, beta_2=0.999),\n",
    "#               loss='sparse_categorical_crossentropy',\n",
    "#               metrics=['accuracy'])\n",
    "# adam_hist = model.fit(train_ds, epochs=5,\n",
    "#           validation_data=test_ds)\n",
    "\n",
    "# losses['Adam'] = adam_hist.history[\"val_loss\"]\n",
    "# make_loss_plots(losses)\n",
    "\n",
    "# SGD with learning rate decay\n",
    "model = build_mlp()\n",
    "model.compile(optimizer=tf.keras.optimizers.SGD(lr=0.5, momentum=0.0, decay=0.01),\n",
    "              loss='sparse_categorical_crossentropy',\n",
    "              metrics=['accuracy'])\n",
    "sgd_decay_hist = model.fit(train_ds, epochs=5,\n",
    "          validation_data=test_ds)\n",
    "\n",
    "losses['SGD with decay'] = sgd_decay_hist.history[\"val_loss\"]\n",
    "make_loss_plots(losses)\n",
    "\n",
    "# Momentum with learning rate decay\n",
    "model = build_mlp()\n",
    "model.compile(optimizer=tf.keras.optimizers.SGD(lr=0.5, momentum=0.9, decay=0.01),\n",
    "              loss='sparse_categorical_crossentropy',\n",
    "              metrics=['accuracy'])\n",
    "sgd_decay_hist = model.fit(train_ds, epochs=5,\n",
    "          validation_data=test_ds)\n",
    "\n",
    "losses['Momentum with decay'] = sgd_decay_hist.history[\"val_loss\"]\n",
    "make_loss_plots(losses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "V3gMgKgltU23"
   },
   "source": [
    "**Exercise:** Play with the parameters of the various optimisers and see how that affects this particular problem.\n",
    "\n",
    "**Note:** There is an element of randomness each time we train our models; this means that we should be careful when making judgements about which methods are best for this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "S5aoZHykLYgW"
   },
   "source": [
    "## Optional extra reading: second-order methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ezeyeKFOOUr1"
   },
   "source": [
    "All of the optimization methods we've looked at in this practical are what we call first order methods: they calculate a straight line estimate of the gradient and take a step in that direction. However, this straight line estimate throws away some useful information about the curvature of our loss surface. In other words, a first order method gives us information about whether or not our loss is increasing or deacresing in a given direction, but a second order method tells us how much it is increasing or decreasing. With this extra information, second order methods can find a minima in fewer steps, and have less trouble with getting stuck in saddle points.\n",
    "\n",
    "Unfortunately, the benefits of second-order methods come at a cost. As you might have guessed from the names, a first-order method calculates the first derivative of our loss function, while a second-order method also calculates the second derivative. Calculating these second derivatives (the [Hessian matrix](https://en.wikipedia.org/wiki/Hessian_matrix)) is computationally expensive, which means that second-order methods are often slower in practice. For this reason, we usually do not see second-order methods in deep learning. However, fast approximations for second-order methods are an active area of research!\n",
    "\n",
    "Examples of second-order methods include [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization) and [BFGS](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VNGBSO8ya2H-"
   },
   "source": [
    "## Conclusion:\n",
    "\n",
    "### What are optimisers, and how do we use them?\n",
    "\n",
    "Optimisers are algorithms that try to find the best item when compared to other items using a given method of comparison (or metric). In deep learning, we use them to try to find the best values for our weights and biases, using a loss function as our metric.\n",
    "\n",
    "### What method should you use?\n",
    "\n",
    "There are **no hard rules** for which methods you should use. It will always depend on your particular model and dataset. However, some guidelines are worth keeping in mind:\n",
    "\n",
    "1. Adam typically works very well in a large number of settings and is usually a good first choice.\n",
    "2. RMSProp can often outperform Adam for RNNs as well as in RL. If you are working in either of these domains, then it might be worth trying RMSProp.\n",
    "3. SGD and SGD with momentum often work just as well as more sophisticated methods like Adam. Don't think that they aren't worth trying out just because they are simple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "f7UcgvT_y5_C"
   },
   "source": [
    "## Tasks\n",
    "\n",
    "1. **[All]** Combine the implementations for momentum and RMSProp to implement Adam. Experiment with the parameters of Adam and compare it to the previous methods we have looked at.\n",
    "2. **[All]** Augment any of the optimization methods with learning rate decay. You can choose any of the learning rate decay methods. Play around with the parameters. Compare the performance of the optimization method with and without learning rate decay.\n",
    "3. **[Optional, Intermediate]** Implement the other two learning rate decay methods and compare all three decay methods with one another.\n",
    "4. **[Optional, Advanced]** Implement Nesterov Momentum. You can read more about it [here](http://ruder.io/optimizing-gradient-descent/).\n",
    "5. **[Optional, Advanced]** Implement Nadam. You can read more about it [here](http://ruder.io/optimizing-gradient-descent/) (it also works particularly well in deep RL).\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "TPiFyajJtDmj"
   },
   "source": [
    "**IMPORTANT: Please fill out the exit ticket form before you leave the practical: https://forms.gle/J4i5wehZPUdggCc29**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1A1GmqWZL_jW"
   },
   "source": [
    "## Extra resources\n",
    "\n",
    "* A [blog post](http://fa.bianp.net/teaching/2018/eecs227at/gradient_descent.html) from Fabian Pedregosa on gradient descent [**Highly Recommended**].\n",
    "* [Distil.pub post](https://distill.pub/2017/momentum/) by Gabriel Goh on why momentum works [**Highly Recommended**].\n",
    "* [Sebastian Ruder's blog](http://ruder.io/optimizing-gradient-descent/) on gradient descent algorithms.\n",
    "* Deep Dive into Deep Learning chapter on [Optimization Algorithms](http://d2l.ai/chapter_optimization/index.html).\n",
    "* Keras optimizer [docs](https://keras.io/optimizers/).\n"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "Khipu Practical 1b: Optimisation for Deep Learning",
   "private_outputs": true,
   "provenance": []
  },
  "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.7.5rc1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}