lecture_notes_0.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "import requests\n",
    "import json\n",
    "from itertools import count\n",
    "from typing import List, NamedTuple, Union\n",
    "from ipywidgets import interact, interactive, fixed, IntSlider\n",
    "from IPython.display import display, HTML\n",
    "from tqdm import tqdm\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.random.seed(20080524)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Multi-layer Perceptrons\n",
    "### University of Cambridge\n",
    "### Engineering Part IIB/Module 4F12: Computer Vision\n",
    "### Deep Learning for Computer Vision\n",
    "\n",
    "*Lecturer: Matthew Johnson*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complex, interconnected networks of neurons, like the human brain, are the only known way of achieving general intelligence. Because of this, many efforts in the area of Artificial Intelligence attempt to recreate the different structures of the brain. Systems called _deep neural nets_ (DNNs) have proven helpful for more targeted AI tasks, like image understanding:\n",
    "\n",
    "[Microsoft Dense Captions API](t{https://aka.ms/densecaptio)\n",
    "\n",
    "Let's see what kind of things a Computer Vision API can tell us about images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyze_image(image_url, api_key):\n",
    "    url = \"https://lecturescv.cognitiveservices.azure.com/vision/v2.1/analyze?visualFeatures=Tags,Description\"\n",
    "    headers = {\"Content-Type\": \"application/json\", \"Ocp-Apim-Subscription-Key\": api_key}\n",
    "    payload = \"{{'url': '{}'}}\".format(image_url)\n",
    "    r = requests.post(url, headers=headers, data=payload)\n",
    "    print(r.status_code, r.reason)\n",
    "    result = json.loads(r.text)\n",
    "\n",
    "    captions = [[\"Caption\", caption[\"text\"], caption[\"confidence\"]] \n",
    "                for caption in result[\"description\"][\"captions\"]]\n",
    "    tags = [[\"Tag\", tag[\"name\"], tag[\"confidence\"]] for tag in result[\"tags\"]]\n",
    "\n",
    "    rows = \"\".join([\"<tr><td>{}</td></tr>\".format(\"</td><td>\".join([str(x) for x in row])) for row in captions + tags])\n",
    "    table = \"<table><tr><th>Property</th><th>Text</th><th>Confidence</th></tr>{}</table>\".format(rows)       \n",
    "\n",
    "    format_string = \"<table><tr><td><div class='img-overlay-wrap'><img src='{}' width='500'/></div></td><td>{}</td></table>\"\n",
    "    html = format_string.format(image_url, table)\n",
    "    \n",
    "    return html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#im_url = \"https://lectures.blob.core.windows.net/media/caption_test_0.jpg\"\n",
    "#im_url = \"https://lectures.blob.core.windows.net/media/caption_test_1.jpg\"\n",
    "im_url = \"https://lectures.blob.core.windows.net/media/caption_test_2.jpg\"\n",
    "with open(\"api_key.txt\") as f:\n",
    "    key = f.read()\n",
    "\n",
    "display(HTML(analyze_image(im_url, key)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DNNs can seem intelligent because they do things we usually associate with human intelligence, like scene understanding and object recognition. In reality, they are highly specialised, trained from vast corpora of data, and constructed from elements that mimic biological neural networks in useful yet incomplete ways. We will explore the underlying mathematics behind DNNs, analyse their structures and design principles, and study state-of-the-art architectures and techniques.\n",
    "\n",
    "Before we dive into the maths, let us perform an (extremely brief) overview of neurons and how they work. Below is a diagram of a synapse:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neurons\n",
    "\n",
    "![Synapse Schematic](https://lectures.blob.core.windows.net/media/synapse.png \"Synapse\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ways in which neurons work together involve complex interactions between electrical signals and neurochemistry that are very much outside the scope of this lecture. That said, for our purposes it is important to understand three key concepts:\n",
    "\n",
    "1. A neuron has both dendrites (inputs) and axons (outputs) which connect it to other neurons\n",
    "2. A neuron will produce an exciting or inhibiting signal along its axons based upon activations from its dendrites in a non-linear way\n",
    "3. Each dendrite contributes to whether a signal is produced in different proportions, which change over time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perceptrons\n",
    "<img style='margin: 10px;width: 30%;' src='https://lectures.blob.core.windows.net/media/rosenblatt.jpg' alt='Frank Rosenblatt' align='right'/>\n",
    "\n",
    "The perceptron algorithm was invented in 1957 at the Cornell Aeronautical Laboratory by Frank Rosenblatt. It was designed to be implemented in hardware for the purpose of image recognition, and was the marvel of the age. The\n",
    "New York Times reported the perceptron to be\n",
    "\n",
    "> the embryo of an electronic computer that [the Navy] expects will be able to walk, talk,\n",
    "> see, write, reproduce itself and be conscious of its existence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roseblatt based his designs upon earlier theoretical work by Warren McCullock and Walter Pitts in 1943 on the Threshold Logic Unit, arguably the first artificial neuron. The idea in both cases was to create a simple mathematical model of a neuron. In the case of the perceptron, this model was then expressed in hardware. The entire machine represented a single neuron in this sense, with 400 inputs from a 20x20 image sensor and a binary output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The math behind perceptrons is fairly straightforward. Given an input $\\mathbf{x}$, the output of a perceptron is calculated as:\n",
    "\n",
    "$$y(\\mathbf{x}) = f(\\mathbf{w}^{\\mathrm{T}}\\mathbf{x})$$\n",
    "\n",
    "where $\\mathbf{w}$ are the weights on the input vector $\\mathbf{x}$ and $f$ is a step function of the form:\n",
    "\n",
    "$$\n",
    "f(x_i) = \\begin{cases}\n",
    "         +1 & x_i \\geq 0 \\\\\n",
    "         -1 & x_i < 0\n",
    "       \\end{cases} \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style='margin: 10px;width:30%' src='https://lectures.blob.core.windows.net/media/perceptron.jpg' align='left' alt='The Mark 1 Perceptron' />\n",
    "\n",
    "The fundamental problem at hand is how to learn the values of $\\mathbf{w}$. For this, we use the _perceptron criterion_:\n",
    "\n",
    "$$E_{\\mathrm{P}}(\\mathbf{w}) = -\\sum_{n \\in \\mathcal{M}}\\mathbf{w}^{\\mathrm{T}}\\mathbf{x}_{n}t_{n}$$\n",
    "\n",
    "where $t_{n} \\in \\{-1, +1\\}$ for negative and positive examples and $\\mathcal{M}$ is the set of misclassified examples.\n",
    "We can use the gradient of this function to update the weights on a per-example basis:\n",
    "\n",
    "$$\\mathbf{w}^{\\tau + 1}=\\mathbf{w}^{\\tau} - \\eta \\nabla E_{\\mathrm{P}}(\\mathbf{w}) = \\mathbf{w}^{\\tau} + \\eta \\mathbf{x}_{n} t_{n}$$\n",
    "\n",
    "where $\\eta$ is the learning rate parameter and $\\tau$ is an integer that indexes the steps of the algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what this looks like in practice.  First, let's create a dataset:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets import BinaryDataset\n",
    "\n",
    "binary_dataset = BinaryDataset.generate_from_normals(200)\n",
    "plt.figure(figsize=(6,6))\n",
    "ax = plt.subplot(111)\n",
    "binary_dataset.plot(ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a 2D dataset generated from two normal distributions. Notice that I have provided a third column which is equal to 1 for each point. This is because we are going to define our decision boundary using the general equation for a line:\n",
    "\n",
    "$$Ax + By + C = 0$$\n",
    "\n",
    "The column of 1s allows $C$ to be learned. This practice is called including a _bias_. Now that we have some points, we want to learn a classifier that distinguishes red from blue. First, let's try creating a classifier by hand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import perceptron\n",
    "\n",
    "def manual_classifier(A, B, C):\n",
    "    plt.figure(figsize=(6,6))\n",
    "    ax = plt.subplot(111)\n",
    "    perceptron.plot_2d_decision_boundary(binary_dataset, perceptron.Snapshot(0, [A,B,C], 0, None), ax, \"Manual Classifier\")\n",
    "\n",
    "interact(manual_classifier, A=(-3.0, 3.0), B=(-3.0, 3.0), C=(-3.0, 3.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What so impressed people about the Perceptron was that it was able to learn the weights automatically just by being shown the data points one by one along with feedback that indicated whether its prediction was right or wrong. Let's see how we would go about implementing the perceptron algorithm we laid out above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_perceptron(dataset: BinaryDataset, num_iterations: int) -> np.ndarray:\n",
    "    num_dims = dataset.num_dims\n",
    "    # initialise the weights\n",
    "    weights = np.zeros((num_dims, 1))  \n",
    "\n",
    "    pos_count = dataset.positive.shape[0]\n",
    "    neg_count = dataset.negative.shape[0]\n",
    "\n",
    "    for i in range(num_iterations):\n",
    "        # select a positive and a negative example\n",
    "        pos = dataset.positive[i % pos_count]\n",
    "        neg = dataset.negative[i % neg_count]\n",
    "\n",
    "        # present the positive example\n",
    "        pos_out = np.dot(pos, weights)\n",
    "\n",
    "        if pos_out < 0:\n",
    "            # if there was a mistake, update the weights\n",
    "            weights = weights + pos.reshape(weights.shape)\n",
    "        \n",
    "        # present the negative example\n",
    "        neg_out = np.dot(neg, weights)\n",
    "        if neg_out >= 0:\n",
    "            # if there was a mistake, update the weights\n",
    "            weights = weights - neg.reshape(weights.shape)\n",
    "\n",
    "    return weights\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have an implementation of the Perceptron algorithm, let's try running it and seeing what happens:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training a perceptron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def training_slideshow_points(snapshots, step):\n",
    "    plt.figure(figsize=(6,6))\n",
    "    ax = plt.subplot(111)\n",
    "    perceptron.plot_2d_decision_boundary(binary_dataset, snapshots[step], ax, \"Perceptron\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snapshots_binary = perceptron.train(binary_dataset, 200)\n",
    "\n",
    "interact(training_slideshow_points, snapshots=fixed(snapshots_binary), step=IntSlider(value=0, min=0, max=len(snapshots_binary)-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the perceptron is able to adjust the decision boundary until it gets all of the points right. The original perceptron was able to operate on 20x20 images and recognise them. To show the perceptron working with this far more difficult task, we're going to use a dataset called [MNIST](http://yann.lecun.com/exdb/mnist/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MNIST\n",
    "\n",
    "MNIST consists of tens of thousands of 28x28 grayscale hand-written digits. \n",
    "\n",
    "![MNIST Examples](https://lectures.blob.core.windows.net/media/mnist_tiled.png \"MNIST\")\n",
    "\n",
    "The images above were gathered from employees of the US Census Bureau and Secondary School students, consisting of\n",
    "approximately 250 individuals. There are 60,000 training images and 10,000 test images (sampled in a disjoint manner\n",
    "from the two sample groups). The images themselves are 28 X 28 unregistered grayscale images. Research has been\n",
    "performed on the dataset since the 1980s and has continued to the present day where it is often used in tutorials for neural net software frameworks. Let's see how our perceptron handles some binary classification challenges from MNIST."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perceptron: 1 vs 0\n",
    "\n",
    "We will begin by looking at how a perceptron predicts whether an image is 1 or 0. Because the weight vector $w$ of a perceptron is the same size as its input, we can visualize these weights as an image. Note how the weights are being set by subtracting or adding images, so that we can see afterimages of specific 1 and 0 images in the early weights. Note also how the weight magnitude increases over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_1_0 = BinaryDataset.mnist(1, 0)\n",
    "snapshots_1_0 = perceptron.train(dataset_1_0, 1000)\n",
    "\n",
    "def training_slideshow_mnist(snapshots, step: int):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    fig.clf()\n",
    "    perceptron.plot_mnist(fig, snapshots[step], snapshots)\n",
    "    plt.show()\n",
    "\n",
    "interact(training_slideshow_mnist, snapshots=fixed(snapshots_1_0), step=IntSlider(value=1, min=1, max=len(snapshots_1_0) - 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perceptron: 2 vs 5\n",
    "\n",
    "Let's look at a more difficult example and see how the perceptron performs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_2_5 = BinaryDataset.mnist(2, 5)\n",
    "snapshots_2_5 = perceptron.train(dataset_2_5, 1000)\n",
    "interact(training_slideshow_mnist, snapshots=fixed(snapshots_2_5), step=IntSlider(value=0, min=1, max=len(snapshots_2_5) - 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining Perceptrons\n",
    "\n",
    "As originally proposed, perceptrons can only answer yes/no questions. However, we can combine perceptrons to tackle multi-class problems using _1 versus all_. Instead of a single vector, each perceptron is now a row in a matrix:\n",
    "\n",
    "$$\\mathbf{y} = f(\\mathbf{W}\\mathbf{x})$$\n",
    "\n",
    "where $\\mathbf{W}$ are the weights on the input vector $\\mathbf{x}$ and $f$ is the same step function, applied on a per-index basis on the output vector:\n",
    "\n",
    "$$\n",
    "f(x_i) = \\begin{cases}\n",
    "         +1 & x_i \\geq 0 \\\\\n",
    "         -1 & x_i < 0\n",
    "       \\end{cases} \\\\\n",
    "$$\n",
    "\n",
    "Each perceptron focuses on a single question: is this input a member of my class, or not? We can then look at multiple outputs and choose the class that corresponds to the classifier that has a positive result. The maths stay much the same, though we are going to change the formulation slightly. Instead of including the bias as an extra dimension in the data, we are going to use a common convention in modern neural net architectures and separate it out as a separate _bias_ term, $\\mathbf{b}$, as follows:\n",
    "\n",
    "$$\\mathbf{y} = f(\\mathbf{W}\\mathbf{x} + \\mathbf{b})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Activation Functions\n",
    "\n",
    "There is a problem as curently stated with our activation function $f$. To see this, look at the result of computing $y$ given the following inputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[0.3, -0.5],[-.4, 1]], dtype='float32')\n",
    "W = np.array([[0.55, -0.12], [-0.32, 0.5]], dtype='float32')\n",
    "b = np.array([0.4, -0.2], dtype='float32')\n",
    "print(\"x=\", x)\n",
    "h = np.dot(x, W.T) + b\n",
    "print(\"h=\", h)\n",
    "y = np.sign(h)\n",
    "print(\"y=\", y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The step function used by the perceptron is just one of a class of functions called _activation functions_. There are many others. One that is a close approximation is the hyperbolic tangent:\n",
    "\n",
    "$$\n",
    "f(x_i) = \\tanh(x_i)\n",
    "$$\n",
    "\n",
    "Unlike the step function, the hyperbolic tangent is fully differentiable, which will come in helpful later. Just as a reminder, that looks like this:\n",
    "\n",
    "$$\n",
    "f'(x_i) = 1 - \\tanh^2(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from activations import tanh, perceptron, compare\n",
    "\n",
    "compare(tanh, perceptron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiperceptron\n",
    "\n",
    "We now have a new formulation for the perceptron:\n",
    "\n",
    "$$\\mathbf{y} = \\tanh(\\mathbf{W}\\mathbf{x} + \\mathbf{b})$$\n",
    "\n",
    "This new formulation, which we will call a _multi-perceptron_, has several advantages. The hyperbolic tangent, in addition to being roughly linear from -1 to 1, saturates to 1 or -1 as well, thus providing a useful _non-linearity_, which was the driving force behind the use of the step function in the original perceptron. It also creates a function that is fully differentiable, which means we can learn the values of $\\mathbf{W}$ via _gradient descent._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent\n",
    "\n",
    "As you may remember, gradient descent is the practice of using the gradient of a function to find the values that correspond to its minimum. Take for example the parabola defined by:\n",
    "\n",
    "$$y(x) = x^2 - 5x + 2$$\n",
    "\n",
    "This can be viewed as an _objective function_ of $x$. We can use its gradient:\n",
    "\n",
    "$$y'(x) = 2x - 5$$\n",
    "\n",
    "to find the value of $x$ that minimises $y(x)$. A step size or _learning rate_ $\\eta$ is used to control how far we move in the direction of the gradient. Let's explore this in practice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gradient_descent\n",
    "\n",
    "t = -2\n",
    "eta = .1\n",
    "steps = 30\n",
    "slideshow = gradient_descent.optimize(t, eta, steps)\n",
    "\n",
    "def gd_slide(step):\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    ax = plt.subplot(111)\n",
    "    gradient_descent.plot_step(ax, slideshow[step])\n",
    "    \n",
    "\n",
    "interact(gd_slide, step=IntSlider(value=0, min=0, max=steps-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Backpropagation\n",
    "\n",
    "The use of gradient descent to discover the parameters of a network of perceptrons is called \\emph{backpropagation}, and was discovered by Paul J. Werbos, who described it in his 1974 thesis _Beyond regression: New tools for prediction and analysis in the behavioural sciences_. In order to use backpropagation here we need two things: a (differentiable) way of producing a desired output value, and an (also differentiable) objective function for scoring that value.\n",
    "\n",
    "The task we are currently interested in performing is _classification_, where we label an input as being a member of one of several classes $C$. One way of achieving this is by generating a probability distribution $P(c~|~\\mathbf{x})$, which defines a distribution over $C$ such that $\\sum_{c \\in C} P(c~|~\\mathbf{x}) = 1$ and where the magnitude of $P(c~|~\\mathbf{x})$ corresponds to the classifier's confidence that $\\mathbf{x}$ belongs to class $c$. A typical choice to produce $P(c~|~\\mathbf{x})$ for a multiperceptron is the Softmax function $s$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    P(c=i~|~\\mathbf{x}) & = s(y_i)                                                \\\\\n",
    "    s(y_i)              & = \\frac{e^{y_i}}{\\sum_j e^{y_j}}                        \\\\\n",
    "    \\mathbf{y}          & = \\tanh(\\mathbf{W}\\mathbf{x} + \\mathbf{b})\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Entropy\n",
    "\n",
    "The final thing we need to train this multiperceptron is an objective function. For a probability distribution, we can use the widely employed cross entropy:\n",
    "\n",
    "$$H(P,Q) = -\\sum_i P(i) \\log Q(i)$$\n",
    "\n",
    "Where $P$ is the target distribution and $Q$ is the current output of the softmax.\n",
    "\n",
    "The cross entropy consists of two components (combined above):\n",
    "    \n",
    "$$H(P,Q) = H(P) + D_{\\mathrm{KL}}(P||Q)$$\n",
    "\n",
    "for discrete distributions $P$ and $Q$, where $H$ is the Shannon entropy:\n",
    "\n",
    "$$H(P) = -\\sum_i P(i)\\log P(i)$$\n",
    "\n",
    "Entropy measures the _information content_ of a distribution. This can be thought of as the expectation of the information content given by sampling from the distribution. A perfectly uniform distribution, with equal probability for all events, has maximum entropy. Think of flipping a coin: if it is fair ($P(Heads) = P(Tails) = 0.5$) then we have no way of predicting what the result will be when it is flipped. However, if the coin were weighted, such that $P(Heads) >> P(Tails)$, then flipping it does not give as much information; we knew it was probably going to be heads.\n",
    "\n",
    "$D_{\\mathrm{KL}}$ is the Kullback-Leibler divergence, which is a common measure of how much one distribution diverges from another:\n",
    "\n",
    "$$D(P||Q) = \\sum_i P(i) \\log \\frac{P(i)}{Q(i)}$$\n",
    "\n",
    "The intuition here is that if $P$ and $Q$ are the same, then the term $\\log \\frac{P(i)}{Q(i)}$ will be zero and thus the divergence is also zero.  We combine them as follows:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "H(P,Q) &= -\\sum_i P(i)\\log P(i) + \\sum_i P(i) \\log \\frac{P(i)}{Q(i)} \\\\\n",
    "&= -\\sum_i P(i)\\log P(i) + \\sum_i P(i)\\log P(i) - \\sum_i P(i)\\log Q(i) \\\\\n",
    "&= -\\sum_i P(i)\\log Q(i)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Let's see how the cross entropy behaves when it compares different probability distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_dists(p):\n",
    "    bins = np.arange(p.shape[1])\n",
    "    fig = plt.figure(figsize=(10, 4))\n",
    "    ax = fig.add_subplot(1, 2, 1)\n",
    "    plt.ylim(0, p.max())\n",
    "    plt.bar(bins, p[0], color='b')\n",
    "    ax.set_xticks([])\n",
    "    ax = fig.add_subplot(1, 2, 2)\n",
    "    plt.ylim(0, p.max())\n",
    "    plt.bar(bins, p[1], color='r')\n",
    "    ax.set_xticks([])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_entropy(p, q):\n",
    "    plogq = p*np.log(q)\n",
    "    return -plogq.sum()\n",
    "dists = np.random.uniform(0, 1, (2, 5)).astype('float32')\n",
    "dists = dists / dists.sum(axis=1, keepdims=True)\n",
    "show_dists(dists)\n",
    "plt.show()\n",
    "print(\"H(p,p) =\", cross_entropy(dists[0],dists[0]))\n",
    "print(\"H(p,q) =\", cross_entropy(dists[0],dists[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimisation\n",
    "\n",
    "We can now write the objective function for a multiperceptron:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    L(D, \\mathbf{W}, \\mathbf{b}) & = \\frac{1}{|D|}\\sum_{d \\in D}H(\\mathbf{t}_d, P(c~|~\\mathbf{x}_d)) \\\\\n",
    "    P(c~|~\\mathbf{x}_d)          & = s(\\mathbf{W}\\mathbf{x}_d + \\mathbf{b})                          \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $D$ is our dataset and $\\mathbf{t}_i$ is a _one-hot_ vector corresponding to the correct label for that data point. In words, we want to minimise the expected cross entropy between the output and the target distribution. Now that we have an objective, we can compute gradients for $\\mathbf{W}$ and $\\mathbf{b}$ and optimize them directly.\n",
    "\n",
    "At first glance the task of computing partial derivatives of $L$ with respect to $\\mathbf{W}$ and $\\mathbf{b}$ seems daunting. However, we have chosen our functions very carefully. Using the chain rule, we first look at partial derivatives with respect to the softmax function, but immediately run into a problem, in that $y_i$ is in both the numerator and denominator of $s$. As such we need to index the derivative:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial s_i}{\\partial y_i} & = s_i(1 - s_i)       \\\\\n",
    "    \\frac{\\partial s_j}{\\partial y_i} & = -s_i s_j,~j \\neq i\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Breaking things out in this way, we can now look at our first partial derivatives for $L$ (for a single data point):\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    L                               & = H(\\mathbf{t}, P(c~|~\\mathbf{x}))                   \\\\\n",
    "    L                               & = -\\sum_{i}t_{i}\\log s_i                             \\\\\n",
    "    \\frac{\\partial L}{\\partial y_i} & = -\\sum_k t_k \\frac{\\partial \\log s_k}{\\partial y_k} \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "We change to $k$ here because $y_i$ shows up in both the numerator and denominator of softmax.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial L}{\\partial y_i} & = -\\sum_k t_k \\frac{1}{s_k} \\frac{\\partial s_k}{\\partial y_k} \\\\\n",
    "                                    & = -t_i(1 - s_i) - \\sum_{k \\neq i} t_k \\frac{1}{s_k}(-s_k s_i)\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "via the chain rule and the derivative of softmax.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial L}{\\partial y_i} & = -t_i + t_i s_i + \\sum_{k \\neq i} t_k s_i \\\\\n",
    "                                    & = s_i\\left(\\sum_k t_k\\right) - t_i         \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Remember that $\\mathbf{t}$ is a one-hot vector, and so we get the final result:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\frac{\\partial L}{\\partial y_i} = s_i - t_i\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "Continuing with the chain rule, we can now look at the partial derivatives for the hyperbolic tangent:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    y_i                               & = \\tanh(h_i)                          \\\\\n",
    "    h_i                               & = \\mathbf{W}\\mathbf{x}_i + \\mathbf{b} \\\\\n",
    "    \\frac{\\partial y_i}{\\partial h_i} & = 1 - y_i^2\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "As we are performing matrix calculus, is it helpful to think about vectors of partial derivatives moving forward, which we will denote $\\delta = \\left\\{\\frac{\\partial L}{\\partial h_i},~\\forall i\\right\\}$. Using $\\delta$ we can compute the partial derivatives for the weights and biases:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial L}{\\partial \\mathbf{W}} & = \\mathbf{x}\\delta^{\\mathrm{T}} \\\\\n",
    "    \\frac{\\partial L}{\\partial \\mathbf{b}} & = \\sum_i \\delta_i\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "By subtracting these values from the current values of $\\mathbf{W}$ and $\\mathbf{b}$, we can perform gradient descent to update the weights and biases and minimise our objective function, with the aim that $P(c~|~x_d) \\approx \\mathbf{t}_d$ for all $d \\in D$.\n",
    "\n",
    "We have everything we need now to code up a multiperceptron:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from activations import dtanh\n",
    "\n",
    "def softmax(x):\n",
    "    x_max = x.max(axis=1,keepdims=True)\n",
    "    x_exp = np.exp(x - x_max) # this avoids numerical issues\n",
    "    x_exp_sum = x_exp.sum(axis=1, keepdims=True)\n",
    "    return x_exp / x_exp_sum\n",
    "\n",
    "def cross_entropy_loss(p, t):\n",
    "    p_of_t = p[np.arange(len(t)), t]\n",
    "    log_prob = np.log(p_of_t)\n",
    "    return -log_prob.mean()\n",
    "\n",
    "def dcross_entropy_loss(p, t):\n",
    "    labels = np.zeros_like(p)\n",
    "    labels[np.arange(len(t)), t] = 1\n",
    "    return (p - labels) / len(t)\n",
    "\n",
    "def linear(x, W, b):\n",
    "    return np.dot(x, W.T) + b\n",
    "\n",
    "def dlinear(x, W, b, dy):\n",
    "    dx = np.dot(dy, W)\n",
    "    dW = np.dot(dy.T, x)\n",
    "    db = dy.sum(axis=0)\n",
    "    return dx, dW, db\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these functions we can compute the derivatives for the parameters of all the perceptrons in one chain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[0.3, -0.5],[-.4, 1]], dtype='float32')\n",
    "W = np.array([[0.55, -0.12], [-0.32, 0.5]], dtype='float32')\n",
    "b = np.array([0.4, -0.2], dtype='float32')\n",
    "\n",
    "print(\"x=\", x)\n",
    "h = linear(x, W, b)\n",
    "print(\"h=\", h)\n",
    "s = tanh(h)\n",
    "print(\"s=\", s)\n",
    "p = softmax(s)\n",
    "print(\"p =\", p)\n",
    "\n",
    "t = np.array([1, 0], dtype='int32')\n",
    "dloss = dcross_entropy_loss(p, t)\n",
    "ds = dtanh(s, dloss)\n",
    "_, dW, db = dlinear(x, W, b, ds)\n",
    "print(\"dloss =\", dloss)\n",
    "print(\"dW =\", dW)\n",
    "print(\"db =\", db)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Batching\n",
    "\n",
    "So far we looked at the equations for only a single data point $\\mathbf{x}$ and its target $\\mathbf{t}$. However, we really want to subtract the value for the full definition of $L$ which incorporates an expectation over all data points:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial L}{\\partial \\mathbf{W}} & = \\frac{1}{|D|}\\sum_{d \\in D} \\mathbf{x}_d \\delta_d^{\\mathrm{T}} \\\\\n",
    "    \\frac{\\partial L}{\\partial \\mathbf{b}} & = \\frac{1}{|D|}\\sum_{d \\in D}\\sum_i \\delta_{di}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "This is a problem, as we need to compute the full gradient before we can update the weights and biases. This is computationally expensive, and so we instead use a technique called _batching_, where we compute the gradient for a subset of the data points, and then update the weights and biases.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\frac{\\partial L}{\\partial \\mathbf{W}} & = \\frac{1}{|D_B|}\\sum_{d \\in D_B} \\mathbf{x}_d \\delta_{d}^{\\mathrm{T}} \\\\\n",
    "    \\frac{\\partial L}{\\partial \\mathbf{b}} & = \\frac{1}{|D_B|}\\sum_{d \\in D_B}\\sum_i \\delta_{di}                    \\\\\n",
    "    D_B                                    & \\subseteq D\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "This is repeated until we have processed all the data points in the dataset. If these subsets are chosen at random, this is called _stochastic gradient descent_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Datasets\n",
    "\n",
    "At this point, we will briefly discuss the correct usage of data. Data should be divided into three sets:\n",
    "\n",
    "1. **Training** - Fitting the model\n",
    "2. **Validation** - Determining hyperparameters\n",
    "3. **Test** - Held out for evaluation.\n",
    "\n",
    "The two commandments of data science can be written in terms of these three sets:\n",
    "\n",
    "1. Thou shalt not train on the test data\n",
    "2. Thou shalt not determine hyperparameters on test data\n",
    "\n",
    "These both derive from the central goal of machine learning: _generalisation_. When we train a model, we want to ensure that it can make reasonable predictions when shown new, unseen data. This is why we take a portion of our labelled data and set it aside for use as test data. If we use the test data to train our model, we have no way of knowing if our model will generalise to that test data, and thus the results are meaningless.\n",
    "\n",
    "The same is true for hyperparameters. If we use the test data to determine hyperparameters, we have no way of knowing if our model will generalise to that test data, and thus the results are meaningless. Thus, we set aside validation data for this purpose. However, we can safely incorporate the validation data when we train the final model as we still have the held-out test data to evaluate the model's performance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training the Multiperceptron\n",
    "\n",
    "Let's see all of this put together into a single algorithm that trains the weights and biases of the multiperceptron:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets import MulticlassDataset\n",
    "\n",
    "def train_multiperceptron(data: MulticlassDataset, learning_rate=0.1, num_epochs=1, batch_size=10, report_frequency=0):\n",
    "    train_data, train_labels = data.train\n",
    "    val_data, val_labels = data.val\n",
    "    num_train = train_data.shape[0]\n",
    "    num_dims = train_data.shape[1]\n",
    "    num_classes = np.max(train_labels) + 1\n",
    "\n",
    "    # Initialize our weights and biases\n",
    "    W = np.zeros((num_classes, num_dims), dtype='float32')\n",
    "    b = np.zeros((num_classes,), dtype='float32')\n",
    "\n",
    "    np.set_printoptions(precision=2, suppress=True)\n",
    "\n",
    "    snapshots: List[Snapshot] = []\n",
    "\n",
    "    def evaluate_model():\n",
    "        # predict the distribution over classes on the validation data\n",
    "        p = forward(val_data, W, b)\n",
    "        \n",
    "        # compute the cross entropy\n",
    "        loss = cross_entropy_loss(p, val_labels)\n",
    "\n",
    "        # get the predicted classes\n",
    "        actual = np.argmax(p, axis=1)\n",
    "\n",
    "        accuracy = (actual == val_labels).sum() / val_data.shape[0]\n",
    "        return accuracy, loss\n",
    "\n",
    "    num_instances = 0\n",
    "    last_accuracy = -1.0\n",
    "    last_instances = 0\n",
    "\n",
    "    for epoch in range(num_epochs):\n",
    "        print(\"epoch\", epoch)\n",
    "        # we permute the training data each epoch so that the minibatches are randomized\n",
    "        perm = np.random.permutation(num_train)\n",
    "\n",
    "        for i in range(0, num_train, batch_size):\n",
    "            # evaluate how we are doing on the validation data\n",
    "            accuracy, loss = evaluate_model()\n",
    "            print(num_instances, ': accuracy =', accuracy, \"loss =\", loss)\n",
    "\n",
    "            # get the minibatch\n",
    "            x = np.asarray(train_data[perm[i:i + batch_size]])\n",
    "            t = np.asarray(train_labels[perm[i:i + batch_size]])\n",
    "\n",
    "            # forward phase\n",
    "            h = linear(x, W, b)\n",
    "            s = tanh(h)\n",
    "            p = forward(x, W, b)\n",
    "            \n",
    "            # compute the loss\n",
    "            loss = cross_entropy_loss(p, t)                \n",
    "\n",
    "            # backpropagation\n",
    "            dloss = dcross_entropy_loss(p, t)\n",
    "            ds = dtanh(s, dloss)\n",
    "            _, dW, db = dlinear(x, W, b, ds)\n",
    "\n",
    "            # update the weights/biases\n",
    "            W -= learning_rate * dW\n",
    "            b -= learning_rate * db\n",
    "            num_instances += x.shape[0]\n",
    "\n",
    "    accuracy, loss = evaluate_model()\n",
    "    print('final: accuracy =', accuracy, \"loss =\", loss)\n",
    "    return W, b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trinary Classification\n",
    "\n",
    "Our multiperceptron is capable of predicting labels for any number of classes, not just the two classes that a generic perceptron can handle. Below we see a dataset containing three classes. We can train our multiperceptron to classify them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets import MulticlassDataset\n",
    "import multiperceptron\n",
    "\n",
    "data_trinary = MulticlassDataset.generate_three_class(200, 100)\n",
    "snapshots_trinary = multiperceptron.train(data_trinary, 0.1, 2, 10)\n",
    "\n",
    "def mp_slide(step):\n",
    "    plt.figure(figsize=(6,6))\n",
    "    ax = plt.subplot(111)\n",
    "    snapshots_trinary[step].plot(ax, data_trinary)\n",
    "\n",
    "interact(mp_slide, step=IntSlider(value=1, min=1, max=len(snapshots_trinary) - 2))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MNIST Classification\n",
    "\n",
    "Now we will look at training a multiperceptron on a harder problem: all ten classes. As was the case with the perceptron we can still view the weights as images, with one weight image per row of the $W$ matrix. As training progresses the weights for each class become positive where that class has pixels, and negative elsewhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_mnist = MulticlassDataset.mnist()\n",
    "snapshots_mnist = multiperceptron.train(dataset_mnist, 0.1, 2, 100)\n",
    "num_snaps = len(snapshots_mnist)\n",
    "\n",
    "def mp_mnist_slide(step):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    multiperceptron.plot_mnist(fig, dataset_mnist, snapshots_mnist, step, num_snaps)\n",
    "\n",
    "interact(mp_mnist_slide, step=IntSlider(value=1, min=1, max=len(snapshots_mnist) - 2))  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uncertainty\n",
    "\n",
    "Aside from the ability to handle multiple classes, the multiperceptron also provides us with a measure of uncertainty. As a point moves closer to or farther away from the decision boundary, the model becomes less or more certain. See what happens as the square moves from the top to the bottom:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_unc = MulticlassDataset.generate_two_class(200, 100)\n",
    "snapshots_unc = multiperceptron.train(dataset_unc, 0.1, 2, 10)\n",
    "\n",
    "W = snapshots_unc[-1].weights\n",
    "b = snapshots_unc[-1].biases\n",
    "start = np.array([0, 3.5], np.float32)\n",
    "end = np.array([2, -3.5], np.float32)\n",
    "dir = end - start\n",
    "\n",
    "plt.rc('font', size=15)\n",
    "def plot_unc(stop):\n",
    "    plt.figure(figsize=(6,6))\n",
    "    ax = plt.subplot(111)\n",
    "    point = start + stop * dir\n",
    "    p = multiperceptron.forward(point.reshape(1, 2), W, b)[0]\n",
    "    multiperceptron.plot_decision_boundaries(ax, dataset_unc.train, W, b, None)\n",
    "    ax.plot(point[0], point[1], 's', c=\"magenta\", ms='15', mew='2.0')\n",
    "    bins = [-4, -3]\n",
    "    ax.bar(bins, p, color=['b', 'r'], width=1)\n",
    "    ax.set_xlim(-5, 5)\n",
    "    ax.set_ylim(-5, 5)\n",
    "    ax.set_title(f\"p = [{p[0]:.2f}, {p[1]:.2f}]\")     \n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "interact(plot_unc, stop=(0.0, 1.0, .05))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Linearity Problem\n",
    "\n",
    "The multiperceptron is a powerful tool, but it has a major flaw: it can only be used on linearly separable data. Below you can see two examples of binary data which is not linearly separable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_xor = MulticlassDataset.generate_xor(200, 100)\n",
    "dataset_bullseye = MulticlassDataset.generate_bullseye(200, 100)\n",
    "snapshots_xor = multiperceptron.train(dataset_xor, 0.05, 5, 10, 10)\n",
    "snapshots_bullseye = multiperceptron.train(dataset_bullseye, 0.05, 5, 10, 10)\n",
    "\n",
    "def plot_xor(fig: plt.Figure, i: int):\n",
    "    ax = fig.add_subplot(1, 2, 1)\n",
    "    snapshots_xor[i].plot(ax, dataset_xor, False, True)\n",
    "    ax.set_xlim(-10, 10)\n",
    "    ax.set_ylim(-10, 10)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"XOR\")\n",
    "\n",
    "    ax = fig.add_subplot(1, 2, 2)\n",
    "    snapshots_xor[i].plot(ax, dataset_bullseye, False, True)\n",
    "    ax.set_xlim(-10, 10)\n",
    "    ax.set_ylim(-10, 10)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"Bullseye\")\n",
    "\n",
    "plt.rc('font', size=15)\n",
    "\n",
    "def xor_slide(step):\n",
    "    fig = plt.figure(figsize=(12, 6))\n",
    "    plot_xor(fig, step)\n",
    "\n",
    "interact(xor_slide, step=IntSlider(value=1, min=1, max=len(snapshots_xor) - 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make no mistake, there is no way to change the multiperceptron to overcome this problem. We will need to change the input data itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Engineering\n",
    "\n",
    "One way to overcome this problem is via _feature engineering_, whereby we engineer a feature transform that takes as input the original data and produces new data which is linearly separable and can be used as input to a multiperceptron. Using some simple features, we can then make these problems solveable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def xor_feature(x: np.ndarray):\n",
    "    return np.stack([x[:, 0], x.prod(1) / np.abs(x[:, 0])], axis=1)\n",
    "\n",
    "def bullseye_feature(x: np.ndarray):\n",
    "    r = np.sqrt((x ** 2).sum(1))\n",
    "    return np.stack([x[:, 0], r], axis=1)\n",
    "\n",
    "dataset_xor_feat = dataset_xor.apply_feature(xor_feature)\n",
    "dataset_bullseye_feat = dataset_bullseye.apply_feature(bullseye_feature)\n",
    "\n",
    "snapshots_xor_feat = multiperceptron.train(dataset_xor_feat, 0.05, 5, 10, 10)\n",
    "snapshots_bullseye_feat = multiperceptron.train(dataset_bullseye_feat, 0.05, 5, 10, 10)\n",
    "\n",
    "def plot_xor_feat(fig: plt.Figure, i: int):\n",
    "    ax = fig.add_subplot(1, 2, 1)\n",
    "    snapshots_xor_feat[i].plot(ax, dataset_xor_feat, False, True)\n",
    "    ax.set_xlim(-10, 10)\n",
    "    ax.set_ylim(-10, 10)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"XOR (x*y)\")\n",
    "\n",
    "    ax = fig.add_subplot(1, 2, 2)\n",
    "    snapshots_bullseye_feat[i].plot(ax, dataset_bullseye_feat, False, True)\n",
    "    ax.set_xlim(-10, 10)\n",
    "    ax.set_ylim(0, 10)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"Bullseye (sqrt(x^2 + y^2))\")\n",
    "\n",
    "plt.rc('font', size=15)\n",
    "\n",
    "def xor_feat_slide(step):\n",
    "    fig = plt.figure(figsize=(12, 6))\n",
    "    plot_xor_feat(fig, step)\n",
    "\n",
    "interact(xor_feat_slide, step=IntSlider(value=1, min=1, max=len(snapshots_xor_feat) - 2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-layer Perceptrons (MLPs)\n",
    "\n",
    "Feature engineering is an incredibly useful tool, but requires human expertise. For low-dimensional datasets like these it can be straightforward, but as data becomes more complex, engineering features becomes increasingly difficult. Ideally, we want to be able to learn the feature transform from the data itself. This is the idea behind _multi-layer perceptrons_.\n",
    "\n",
    "We will add a multiperceptron as an intermediary step between the input data and our classifier. This multiperceptron has a different task: to learn a feature transform that makes the data linearly separable. When used in this way, we will refer to a multiperceptron as a _hidden fully-connected layer_. The output of one or more hidden layers is fed into the classifier, which is now called the _output layer_. The resulting network is called a _multi-layer perceptron_ or MLP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyTorch\n",
    "\n",
    "[PyTorch](https://pytorch.org/) is a framework for use in building gradient-based optimization systems like neural nets. It is a Python-first library, built on top of automatic differentiation technology. What this means in practice is that we can write clean, easy to understand Python code and the framework will automatically figure out how to compute the derivative of our calculation for use in backpropagation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "\n",
    "x = torch.rand(4, 2, dtype=torch.float32)\n",
    "print(x)\n",
    "W = torch.rand(3, 2, dtype=torch.float32, requires_grad=True)\n",
    "print(W)\n",
    "b = torch.zeros(3, dtype=torch.float32, requires_grad=True)\n",
    "print(b)\n",
    "h = F.linear(x, W, b)\n",
    "print(h)\n",
    "s = h.tanh()\n",
    "p = F.softmax(s, dim=1)\n",
    "grad = torch.tensor([[1,0,0],[0,1,0],[0,0,1],[0,1,0]], dtype=torch.float32)\n",
    "p.backward(grad)\n",
    "print(W.grad)\n",
    "print(b.grad)\n",
    "#the gradients are now stored in W.grad, b.grad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While we can do the math ourselves, we can also use existing building blocks to build a neural net module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MultiLayerPerceptron(nn.Module):\n",
    "    def __init__(self, input_dim, hidden_dim, output_dim):\n",
    "        nn.Module.__init__(self)\n",
    "        self.linear0 = nn.Linear(input_dim, hidden_dim)\n",
    "        self.linear1 = nn.Linear(hidden_dim, output_dim)\n",
    "\n",
    "    def __call__(self, x):\n",
    "        h = self.linear0(x)\n",
    "        s = h.tanh()\n",
    "        y = self.linear1(s)\n",
    "        \n",
    "        return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code to train this module is similar to what we've seen previously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_multilayer_perceptron(dataset, net: nn.Module, criterion, optimizer, num_epochs = 1, batch_size = 100):\n",
    "    train_data, train_labels = dataset.train\n",
    "    val_data, val_labels = dataset.val \n",
    "    \n",
    "    train_data = torch.from_numpy(train_data)\n",
    "    train_labels = torch.from_numpy(train_labels)\n",
    "    val_data = torch.from_numpy(val_data)\n",
    "    val_labels = torch.from_numpy(val_labels)\n",
    "    \n",
    "    num_train = train_data.shape[0]\n",
    "    \n",
    "    np.set_printoptions(precision=2, suppress=True)\n",
    "    num_instances = 0    \n",
    "    \n",
    "    def compute_accuracy(output, expected_labels):\n",
    "        _, actual_labels = torch.max(output, -1)\n",
    "        num_correct = (actual_labels == expected_labels).sum().item()\n",
    "        total = len(actual_labels)\n",
    "        return num_correct / total\n",
    "   \n",
    "    def evaluate_model():\n",
    "        # put the model in \"evaluation\" mode to ensure no gradients leak from the validation data\n",
    "        net.eval()\n",
    "        x = val_data\n",
    "        t = val_labels\n",
    "        y = net(x)\n",
    "        loss = criterion(y, t).item()\n",
    "        accuracy = compute_accuracy(y, t)\n",
    "        \n",
    "        # put the model back in \"training\" mode\n",
    "        net.train()\n",
    "        return accuracy, loss\n",
    "      \n",
    "    for epoch in range(num_epochs):\n",
    "        print(\"epoch\", epoch)\n",
    "        \n",
    "        # permute the training indices to ensure random batches\n",
    "        perm = np.random.permutation(num_train)\n",
    "        \n",
    "        for i in range(0, num_train, batch_size):\n",
    "            accuracy, loss = evaluate_model()\n",
    "            print(num_instances, \": accuracy =\", accuracy, \"loss =\", loss)                \n",
    "                \n",
    "            # get the minibatch\n",
    "            x = train_data[perm[i:i + batch_size]]\n",
    "            t = train_labels[perm[i:i + batch_size]]\n",
    "            \n",
    "            # zero out the gradients\n",
    "            optimizer.zero_grad()\n",
    "\n",
    "            # forward pass\n",
    "            y = net(x)\n",
    "            \n",
    "            # compute the loss\n",
    "            loss = criterion(y, t)\n",
    "            accuracy = compute_accuracy(y, t)\n",
    "            \n",
    "            # backpropagation\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "           \n",
    "            num_instances += x.shape[0]\n",
    "                    \n",
    "    accuracy, loss = evaluate_model()\n",
    "    print(\"final: accuracy =\", accuracy, \"loss =\", loss)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how the MLP approaches solving the XOR problem. On the left is the input space, with each pixel colored based upon its assigned class. In the middle we see the three-dimensional latent space. You can rotate it using the \"azimuth\" control, and turn on the decision planes for the two classes. Finally, we see accuracy and loss over the course of training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mlp\n",
    "\n",
    "model_xor = mlp.MultiLayerPerceptron(2, 3, 2)\n",
    "if os.path.exists(\"xor.results\"):\n",
    "    results = torch.load(\"xor.results\")\n",
    "    print(\"xor.results loaded\")\n",
    "else:\n",
    "    criterion = nn.CrossEntropyLoss()\n",
    "    optimizer = torch.optim.Adam(model_xor.parameters(), lr=0.01)\n",
    "    results = mlp.train(dataset_xor, model_xor, criterion, optimizer, num_epochs=5, batch_size=5, report_frequency=10)\n",
    "    torch.save(results, \"xor.results\")\n",
    "\n",
    "snapshots_xor_mlp = results['snapshots']\n",
    "\n",
    "plt.rc('font', size=15)\n",
    "\n",
    "def mlp_xor_slide(step, azimuth, planes):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    model_xor.load_state_dict(snapshots_xor_mlp[step].state_dict)\n",
    "    mlp.plot_xor(fig, snapshots_xor_mlp, step, model_xor, dataset_xor, azimuth=azimuth, planes=planes) \n",
    "\n",
    "interact(mlp_xor_slide, step=IntSlider(value=1, min=1, max=len(snapshots_xor_mlp)-1), azimuth=(0,360), planes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the MLP uses the third dimension to make the classes linearly separable. It uses the same trick to solve the bullseye problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_bullseye = mlp.MultiLayerPerceptron(2, 3, 2)\n",
    "if os.path.exists(\"bullseye.results\"):\n",
    "    results = torch.load(\"bullseye.results\")\n",
    "    print(\"bullseye.results loaded\")\n",
    "else:\n",
    "    criterion = nn.CrossEntropyLoss()\n",
    "    optimizer = torch.optim.Adam(model_bullseye.parameters(), lr=0.01)\n",
    "    results = mlp.train(dataset_bullseye, model_bullseye, criterion, optimizer, num_epochs=5, batch_size=5, report_frequency=10)\n",
    "    torch.save(results, \"bullseye.results\")\n",
    "\n",
    "snapshots_bullseye_mlp = results['snapshots']\n",
    "\n",
    "plt.rc('font', size=15)\n",
    "\n",
    "def mlp_bullseye_slide(step, azimuth, planes):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    model_bullseye.load_state_dict(snapshots_bullseye_mlp[step].state_dict)\n",
    "    mlp.plot_xor(fig, snapshots_bullseye_mlp, step, model_bullseye, dataset_bullseye, azimuth=azimuth, planes=planes) \n",
    "\n",
    "interact(mlp_bullseye_slide, step=IntSlider(value=1, min=1, max=len(snapshots_bullseye_mlp)-1), azimuth=(0,360), planes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Understanding the Latent Space\n",
    "\n",
    "The latent spaces learned by MLPs can be difficult to understand. However, one useful technique is to see how values change when we use the differential nature of the network explore how input changes cause output changes. Instead of computing partial derivatives with respect to the weights and biases of the network, we compute them with respect to the input data itself. By changing the target class and then updating the input data to minimise this new objective (while keeping the model fixed) we can get a sense for how the model makes its decisions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = torch.load(\"xor.results\")\n",
    "snapshots_xor_mlp = results['snapshots']\n",
    "model_xor = mlp.MultiLayerPerceptron(2, 3, 2)\n",
    "model_xor.load_state_dict(snapshots_xor_mlp[-1].state_dict)\n",
    "criterion = nn.CrossEntropyLoss()\n",
    "\n",
    "x = dataset_xor.train.values[:1].copy()\n",
    "t = 1 - dataset_xor.train.labels[:1]\n",
    "snapshots_xor_explore, examples_xor = mlp.explore(x, t, model_xor, criterion)\n",
    "\n",
    "plt.rc(\"font\", size=15)\n",
    "\n",
    "def xor_explore_slide(step, azimuth):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    mlp.plot_xor(fig, snapshots_xor_explore, step, model_xor, dataset_xor, examples_xor, azimuth=azimuth)\n",
    "    \n",
    "interact(xor_explore_slide, step=IntSlider(value=1, min=1, max=len(snapshots_xor_explore)-1), azimuth=(0,360))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MNIST: Multi-layer Perceptron\n",
    "\n",
    "Let's look at training a multi-layer perceptron on the same problem. Here is a diagram of a neural net:\n",
    "\n",
    "![MNIST MLP](https://lectures.blob.core.windows.net/media/mlp.svg)\n",
    "\n",
    "We will be using these network diagrams moving forward, so it is worth spending a bit of time understanding it. The layers of the MLP are represented as vertical rectangles. The circles within them represent relative numbers of perceptrons per layer. Above each rectangle is an identifier (\"F\" here for **F**ully Connected) and a number indicating the number of perceptrons. Below is the dimensionality of the output. Modifications of the output, such as the hyperbolic tangent, are shown as superimposed circles over the arrow connecting two layers.\n",
    "\n",
    "The network show here has a hidden layer dimension of 64, meaning that first the 784-dimensional image is embedded in a 64-dimensional space, and then the multiperceptron determines the classification boundary as a hyperplane in 64 dimensions. We can view the weights in the first layer the same way we have done so far, and can view the weights in the second layer as a $5 \\times 2$ grid of $8x8$ images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_mnist = mlp.MultiLayerPerceptron(784, 64, 10)\n",
    "\n",
    "if os.path.exists(\"mnist.results\"):\n",
    "    results = torch.load(\"mnist.results\")\n",
    "    print(\"mnist.results loaded\")\n",
    "else:\n",
    "    criterion = nn.CrossEntropyLoss()\n",
    "    optimizer = torch.optim.SGD(model_mnist.parameters(), lr=0.1)\n",
    "    results = mlp.train(dataset_mnist, model_mnist, criterion, optimizer, num_epochs=5)\n",
    "    torch.save(results, \"mnist.results\")\n",
    "\n",
    "snapshots_mlp_mnist = results[\"snapshots\"]\n",
    "plt.rc('font', size=15)\n",
    "\n",
    "def mlp_mnist_slide(step):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    mlp.plot_mnist(fig, snapshots_mlp_mnist, step, model_mnist)\n",
    "\n",
    "interact(mlp_mnist_slide, step=IntSlider(value=1, min=1, max=len(snapshots_mlp_mnist)-2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Attacking an MLP Model\n",
    "\n",
    "As we did before, we can see how an input image needs to change to cause the network to predict a different class. However, unlike the previous example where the input and latent spaces were low dimensional, these images are embedded in a 784-dimensional space. If we alter the objective function to additionally minimise the distance in the input space:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    L'(\\mathbf{W}, \\mathbf{b}) = L({d}, \\mathbf{W}, \\mathbf{b}) + \\lambda \\sum_{d \\in D} \\left\\lVert \\mathbf{x}_d - \\mathbf{x}'_d \\right\\rVert^2\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "where $\\mathbf{x}'_d$ is the current altered image, can we turn this 4 into a 9?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = torch.load(\"mnist.results\")\n",
    "print(\"mnist.results loaded\")\n",
    "snapshots_mlp_mnist = results['snapshots']\n",
    "model_mnist = MultiLayerPerceptron(784, 64, 10)\n",
    "model_mnist.load_state_dict(snapshots_mlp_mnist[-1].state_dict)\n",
    "\n",
    "def criterion(x0, x, y, t):\n",
    "    return F.cross_entropy(y, t) + 0.01 * F.mse_loss(x, x0, reduction='sum')\n",
    "\n",
    "x = dataset_mnist.train.values[2:3].copy()\n",
    "t0 = dataset_mnist.train.labels[2:3].copy().item()\n",
    "t1 = dataset_mnist.train.labels[4:5].copy().item()\n",
    "\n",
    "snapshots_mnist_attack, examples_attack, labels_attack = mlp.attack(x, t0, t1, model_mnist, criterion)\n",
    "\n",
    "plt.rc('font', size=15)\n",
    "\n",
    "def mnist_attack_slide(step):\n",
    "    fig = plt.figure(figsize=(12, 4))\n",
    "    mlp.plot_mnist_explore(fig, snapshots_mnist_attack, examples_attack, labels_attack, step)\n",
    "\n",
    "interact(mnist_attack_slide, step=IntSlider(value=1, min=1, max=len(snapshots_mnist_attack)-2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Momentum\n",
    "\n",
    "Before we move on to more complex deep learning structures, we will look at a few tricks of the trade. The first is a concept called _momentum_. A typical update in gradient descent is:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\mathbf{W}^{\\tau} \\leftarrow \\mathbf{W}^{\\tau - 1} - \\eta \\frac{\\partial L}{\\partial \\mathbf{W}^{\\tau - 1}}\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "where $\\eta$ is the learning rate and $\\tau$ is the current training step. During learning it is often the case that there is a clear gradient direction, and we can move quickly along it. Conversely, if there is disagreement between subsequent gradients, we may want to move more slowly. We can achieve this by adding a momentum term:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\nabla \\mathbf{W}^{\\tau} & = \\frac{\\partial L}{\\partial \\mathbf{W}^{\\tau - 1}} + \\epsilon \\nabla \\mathbf{W}^{\\tau - 1} \\\\\n",
    "    \\mathbf{W}^{\\tau}        & \\leftarrow \\mathbf{W}^{\\tau - 1} - \\eta \\nabla \\mathbf{W}^{\\tau}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $\\epsilon$ is the momentum parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = 30\n",
    "slideshow = gradient_descent.optimize(-2, 0.15, steps)\n",
    "slideshow_m = gradient_descent.optimize(7, 0.15, steps, momentum=0.4)\n",
    "\n",
    "def gd_mom_slide(step):\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    ax = plt.subplot(111)\n",
    "    gradient_descent.plot_step(ax, slideshow[step])\n",
    "    gradient_descent.plot_step(ax, slideshow_m[step], \"g^\", \"g\", \"dx+mom\")\n",
    "    ax.legend(loc='upper center')\n",
    "    \n",
    "\n",
    "interact(gd_mom_slide, step=IntSlider(value=0, min=0, max=steps-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Rate Updates\n",
    "\n",
    "As the network converges to a local minimum, it is often necessary to make smaller adjustments in the weight values, which translates to smaller step sizes during weight updates. A simple way of doing this is by using a learning rate that changes over time in a stepped fashion:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\eta^{\\tau} = \\eta^{0}\\gamma^{\\lfloor \\tau/\\sigma \\rfloor }\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "where the initial learning rate $\\eta^{0}$ is multiplied by $\\gamma$ every $\\sigma$ steps.\n",
    "\n",
    "Another method is to continuously alter the learning rate in inverse proportion to the number of steps:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\eta^{\\tau} = \\frac{\\eta^{0}}{(1 + \\gamma\\tau)^{\\rho}}\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "where $\\gamma$ and $\\rho$ control the speed of learning rate decay. Let's see what this looks like over the course of training:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "t = np.arange(0, 2900, 100)\n",
    "eta = 0.05\n",
    "stepped = eta * np.power(0.7, np.floor(t / 1000))\n",
    "inverse = eta / np.power(1 + 0.1*t, 0.1)\n",
    "\n",
    "plt.plot(t, stepped, label=\"Stepped LR\")\n",
    "plt.plot(t, inverse, label=\"Exponential Decay\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another option is the popular Adam algorithm, which uses a running average of the gradient and its square to compute the learning rate:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\mathbf{m}^{\\tau}       & = \\beta_1 \\mathbf{m}^{\\tau - 1} + (1 - \\beta_1) \\frac{\\partial L}{\\partial \\mathbf{W}^{\\tau - 1}}                \\\\\n",
    "    \\mathbf{v}^{\\tau}       & = \\beta_2 \\mathbf{v}^{\\tau - 1} + (1 - \\beta_2) \\left(\\frac{\\partial L}{\\partial \\mathbf{W}^{\\tau - 1}}\\right)^2 \\\\\n",
    "    \\hat{\\mathbf{m}}^{\\tau} & = \\frac{\\mathbf{m}^{\\tau}}{1 - \\beta_1^{\\tau}} ~~~~~~\n",
    "    \\hat{\\mathbf{v}}^{\\tau} = \\frac{\\mathbf{v}^{\\tau}}{1 - \\beta_2^{\\tau}}                                                                     \\\\\n",
    "    \\mathbf{W}^{\\tau}       & \\leftarrow \\mathbf{W}^{\\tau} - \\frac{\\alpha \\hat{\\mathbf{m}}^{\\tau}}{\\sqrt{\\hat{\\mathbf{v}}^{\\tau}} + \\epsilon}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Vanishing Gradients\n",
    "\n",
    "There are limits to the level of transformation that can be performed by a single layer of an MLP. As such, it is tempting to add more layers to the network to increase its expressive power. However, as we add more layers, we run into a problem: the gradients become very small, and thus the weights are not updated. This is called the \\emph{vanishing gradient problem}. Here are the derivatives of the most common activation functions from early neural nets: the hyperbolic tangent and sigmoid functions:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\tanh(x)                             & = \\frac{e^{x} - e^{-x}}{e^{x} + e^{-x}} \\\\\n",
    "    \\frac{\\partial \\tanh(x)}{\\partial x} & = 1 - \\tanh^2(x)                        \\\\\n",
    "    s(x)                                 & = \\frac{1}{1 + e^{-x}}                  \\\\\n",
    "    \\frac{\\partial s(x)}{\\partial x}     & = s(x)(1 - s(x))                        \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The derivative of both is always less than 1, and thus the gradient magnitude will decrease as we add more layers, as can be seen with the gradient for the sigmoid function (note the log scale):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlp.gradients(\"sigmoid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One alternative to these activation functions is the \\emph{rectified linear unit} or ReLU:\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    f(x_i) = \\begin{cases}\n",
    "        x_i & x_i \\geq 0 \\\\\n",
    "        0   & x_i < 0\n",
    "    \\end{cases} \\\\\n",
    "\\end{equation*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from activations import relu\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "x_values = np.linspace(-5, 5, 100)\n",
    "plt.ylim(-5, 5)\n",
    "plt.plot(x_values, relu(x_values))\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aside from the discontinuity at zero it is trivially differentiable and, importantly, does not saturate. This means that the gradient will not vanish as we add more layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlp.gradients(\"relu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reducing Model Capacity\n",
    "\n",
    "At the moment our input transformations are performed by a fully connected layer, where each neuron in the output is connected to every input. While these layers are very expressive and can handle a wide range of potential inputs, they also require many parameters. Let's look at our example network from earlier again:\n",
    "\n",
    "![MNIST MLP](https://lectures.blob.core.windows.net/media/mlp.svg)\n",
    "\n",
    "The first layer (from the 784 pixels to the 64 dimensional latent space) requires 50,240 (!) parameters. As a rule, the more parameters you have, the more capacity your network has to overfit to a dataset through memorisation. This hurts our overarching goal of generalisation. As such, we want to reduce the number of parameters in our network (if possible). One way of doing this is by taking advantage of known structure in the input data.\n",
    "\n",
    "To see what we mean by structure, take a look at the images below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets import load_coco\n",
    "\n",
    "dataset_coco = load_coco(\"minival\")\n",
    "images = dataset_coco[\"images\"][:20]\n",
    "\n",
    "def plot(fig: plt.Figure, frame: int):\n",
    "    image = images[frame]\n",
    "    sample = np.random.uniform(0, 1, size=image.shape)\n",
    "    scramble = image.copy().reshape(-1)\n",
    "    np.random.shuffle(scramble)\n",
    "    scramble = scramble.reshape(image.shape)\n",
    "\n",
    "    ax = fig.add_subplot(131)\n",
    "    ax.imshow(sample, interpolation='nearest', vmin=0, vmax=1)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"Sample\")\n",
    "\n",
    "    ax = fig.add_subplot(132)\n",
    "    ax.imshow(image, interpolation='nearest', vmin=0, vmax=1)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"Structured\")\n",
    "\n",
    "    ax = fig.add_subplot(133)\n",
    "    ax.imshow(scramble, interpolation='nearest', vmin=0, vmax=1)\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_title(\"Unstructured\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "\n",
    "def image_structure_slide(step):\n",
    "    fig = plt.figure(figsize=(9, 4))\n",
    "    plot(fig, step)\n",
    "\n",
    "interact(image_structure_slide, step=(0,len(images)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the left we see a sample from the space of all $256 \\times 256$ RGB images. In the middle, we have a natural image. Finally, on the right we have the pixel data from the structure image randomly shuffled. Clearly, natural images are a significant subset of all possible images (i.e. the samples on the left). However, to the MLP, a structured image looks identical to the unstructured pixels on the right. It has to learn the underlying meaning of images as part of the learning process, and much of the capacity is used for this purpose. If we can find a way to help a neural net \"see\" by exploiting image structure, we can significantly reduce model capacity and thus improve generalisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.11"
  },
  "livereveal": {
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}