02_neural_networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creative machine learning - Neural networks\n",
    "\n",
    "### Author: Philippe Esling (esling@ircam.fr)\n",
    "\n",
    "In this course we will cover\n",
    "1. A [quick introduction](#intro) on the principles of neural networks\n",
    "2. An implementation for a [single neuron](#neuron) in Numpy and JAX.\n",
    "3. An exercise on [multi-layer perceptron (MLP)](#mlp) through manual derivation.\n",
    "4. An introduction on [using Pytorch](#pytorch) for defining networks\n",
    "5. An exercise on [audio classification](#audio) using an MLP with Pytorch."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"intro\"></a>\n",
    "# Introducing neural networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will cover more advanced models known as *neural networks*. The tutorial starts by performing a simple **single neuron** discrimination of two random distributions. We will exhibit the manual implementation using Numpy, and then simplify it with JAX. Then, we will study the typical **XOR problem** by using a more advanced 2-layer **perceptron**. Finally, we generalize the use of neural networks in order to perform classification on a given set of audio files, using the PyTorch library, which will provide simplified implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to use relatively _low-level_ libraries to perform the first exercises (implementing your own neurons). To observe this idea in simple setups, we are going to use the `numpy` library and also initialize the homemade course library `cml` and style for future plotting and exercise. We also set the random generator to a fixed point with `rng = np.random.RandomState(1)`, to ensure reproducibility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Base imports\n",
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from cml.plot import initialize_bokeh\n",
    "from cml.panel import initialize_panel\n",
    "from jupyterthemes.stylefx import set_nb_theme\n",
    "from bokeh.io import show\n",
    "initialize_bokeh()\n",
    "initialize_panel()\n",
    "set_nb_theme(\"onedork\")\n",
    "rng = np.random.RandomState(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, to simplify your work, we provide sets of functions that allows fast problem definition and plotting functionnalities (from our `cml` library).\n",
    "\n",
    "  |**File**|*Explanation*|\n",
    "  |-------:|:---------|\n",
    "  |`scatter_boundary`|Plots the decision boundary of a single neuron with 2-dimensional inputs|\n",
    "  |`scatter_classes`|Plots (bi-dimensionnal) input patterns|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from cml.plot import scatter_classes, scatter_boundary, center_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that a single neuron is only able to learn _linearly separable_ problems. To produce such classes of problems, we provide a script that draw a set of random 2-dimensional points, then choose a random line in this space that will act as the linear frontier between 2 classes (hence defining a linear 2-class problem). The variables that will be used by your code are the following.  \n",
    "\n",
    "```Python\n",
    "y_class       # classes of the observqtions \n",
    "x_inputs      # 2 x n final matrix of random input observations\n",
    "weights       # 2 x 1 vector of neuron weights\n",
    "bias          # 1 x 1 vector of bias\n",
    "```\n",
    "\n",
    "You can execute the code below to see our simple classification problem. (Note that running the same cell multiple times produces a different starting dataset). In order to have a well-defined classification problem, we can rely on the `make_blobs` function provided by `scikit-learn` in the `sklearn.datasets` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.datasets import make_blobs\n",
    "# Properties of the problem\n",
    "n_observations = 200\n",
    "noise = 0.2\n",
    "c1_center = [-2, -1]\n",
    "c2_center = [2, 1]\n",
    "# Create points\n",
    "x_coords, y_class = make_blobs(n_samples=n_observations, centers=[c1_center, c2_center], n_features=2, cluster_std=0.55)\n",
    "x_inputs = x_coords + (noise * np.random.randn(n_observations, 2))\n",
    "# Plot the corresponding pattern\n",
    "plot = center_plot(scatter_classes(x_inputs, y_class))\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous example does not allow us to know the exact values (slope and bias) of our ground truth separation line. Hence, we can define a more complex separation problem (with points laying almost right on the separation fronteer), but with known values for our ground truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Number of points to generate\n",
    "n_observations = 100;\n",
    "# Generate 2-dimensional random points\n",
    "x_inputs = np.random.rand(int(n_observations), 2) * 2 - 1;\n",
    "# Slope of separating line\n",
    "sep_slope = np.log(np.random.rand() * 10);\n",
    "sep_bias = np.random.rand() * 2 - 1;\n",
    "# Create the indexes for a two-class problem\n",
    "y_class = (x_inputs[:, 1] - x_inputs[:, 0] * sep_slope - sep_bias > 0) * 1;\n",
    "# Plot the corresponding pattern\n",
    "plot = center_plot(scatter_classes(x_inputs, y_class))\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"neuron\"></a>\n",
    "## Single neuron\n",
    "\n",
    "For the first parts of the tutorial, we will perform the simplest classification model possible in a neural network setting, a single neuron. We briefly recall here that; given an input vector $ \\mathbf{x} \\in \\mathbb{R}^{n} $, a single neuron computes the function  \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\bar{y}=\\phi\\left(\\sum_{i = 1}^{n}w_{i}.x_{i} + b\\right)\n",
    "\\label{eq1}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "with $ \\mathbf{w} \\in \\mathbb{R}^{n} $ a weight vector, $ b $ a bias and $ \\phi\\left( \\cdot \\right) $ an *activation function*. Therefore, if we consider the *threshold* activation function ($ \\mathbb{I}_{0}\\left(x\\right)=1 $ if $ x \\geq 0$), a single neuron simply performs an *affine transform* and then a *linear* discrimination of the space. \n",
    "\n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #410819; border-color: #cb2e47\">\n",
    "\n",
    "> **Important note**\n",
    "> The following implementation tries to remain as close as possible to the original neuron model by McCulloch & Pitts, notably by using the _threshold activation_ function. Although this leads to a very simple implementation (akin to using an _identity_ activation, note that it is not to be used afterwards.\n",
    "\n",
    "</div>\n",
    "\n",
    "***\n",
    "\n",
    "As we will see later on, a **neural network** will simply be composed of _layers_ of these neurons, which produce successive computations\n",
    "\n",
    "<img src=\"images/02_feedforward_nn.png\" align=\"center\"/>\n",
    "\n",
    "Geometrically, a single neuron computes an hyperplane that separates the space. In order to learn, we have to adjust the weights and know \"how much wrong we are\". To do so, we consider that we know the desired output $ y_{j} $ of a system for a given example $ \\mathbf{x}_{j} $ (eg. a predicted value for a regression system, a class value for a classification system). Therefore, we define the MSE loss function $ \\mathcal{L}_{\\mathcal{D}} $ over a whole dataset as\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathcal{L}_{\\mathcal{D}}=\\sum_{j=1}^{|\\mathcal{D}|}\\left\\Vert \\bar{y}_{j}-y_{j}\\right\\Vert ^{2}\n",
    "\\label{eq2}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "with $|\\mathcal{D}|$ being the size of our dataset. In order to know how to change the weights based on the value of the errors, we need to now \"how to change it to make it better\". Therefore, we should compute the sets of derivatives of the error given each parameter\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\Delta\\bar{\\mathbf{w}}=\\left(\\frac{\\delta\\mathcal{L}_{\\mathcal{D}}}{\\delta w_{1}},\\ldots,\\frac{\\delta\\mathcal{L}_{\\mathcal{D}}}{\\delta w_{n}}\\right)\n",
    "\\label{eq3}\n",
    "\\end{equation}\n",
    "$$ \n",
    "\n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #192841; border-color: #779ecb\">\n",
    "\n",
    "> ### Exercise (course)\n",
    ">   1. Perform the derivatives of the output given a single neuron\n",
    ">   2. Perform the derivatives for the bias as well\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training your own neuron\n",
    "\n",
    "We will start by training a single neuron to learn how to perform this discrimination with a linear problem (so that a single neuron is enough to solve it). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Learning rate\n",
    "eta = 0.005;\n",
    "# Weight decay\n",
    "lambda_r = 0.1\n",
    "# Number of epochs\n",
    "n_epochs = 50\n",
    "# Initialize the weights\n",
    "weights = np.random.randn(1, 2);\n",
    "bias = np.random.randn(1, 1);\n",
    "# Save the weight history for plotting\n",
    "weights_history = np.zeros((n_epochs + 1, 2));\n",
    "bias_history = np.zeros((n_epochs + 1, 1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you need to update the following code loop to ensure that your neuron learns to separate between the classes \n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #192841; border-color: #779ecb\">\n",
    "\n",
    "> ### Exercise (course)\n",
    ">   1. Update the loop so that it computes the forward propagation error\n",
    ">   2. Update the loop to perform learning (based on back-propagation)\n",
    ">   3. Run the learning procedure, which should produce a result similar to that displayed on the website\n",
    ">   4. Perform multiple re-runs by **tweaking the hyperparameters** (learning rate, weight decay)\n",
    ">   5. What observations can you make on the learning process?\n",
    ">   6. (Optional) Change the input patterns, and confirm your observations.\n",
    ">   6. (Optional) Incorporate the bias in the weights to obtain a **vectorized** code.\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "verbose = False\n",
    "weights_history[0] = weights\n",
    "bias_history[0] = bias\n",
    "# Update loop\n",
    "for i in range(n_epochs):\n",
    "    \n",
    "    ######################\n",
    "    # Solution  \n",
    "    y_bar = np.dot(weights, x_inputs.T) + bias > 0\n",
    "    if np.sum(np.abs(y_bar - y_class)) == 0:\n",
    "        break\n",
    "    error = y_class - y_bar\n",
    "    weights = weights + eta * (np.dot(error, x_inputs))\n",
    "    bias = bias + eta * (np.sum(error))\n",
    "    ######################\n",
    "    \n",
    "    weights_history[i + 1] = weights\n",
    "    bias_history[i + 1] = bias\n",
    "    \n",
    "    if (verbose):\n",
    "        print('%2d. error = %f, weights = %f, %f, %f' % (i, np.sum(np.abs(error)) / n_observations, bias[0, 0], weights[0, 0], weights[0, 1]))\n",
    "plot = center_plot(scatter_boundary(x_inputs, y_class, weights_history, bias_history, i + 1))\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neuron implementation in `JAX`\n",
    "\n",
    "As seen in the previous course, the `NumPy` implementation requires us to perform manual differentiation of our loss function to understand how to update the parameters. However, the recent [`JAX`](https://github.com/google/jax) library provides **automatic differentiation** on top of the `NumPy` library. If you have not yet completed the set of [tutorials](https://jax.readthedocs.io/en/latest/), we strongly encourage you to do so. We will introduce a few more functions compared to last time, but we try to keep the implementation as simple as possible on purpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import jax.numpy as jnp\n",
    "from jax import grad, jit, random\n",
    "key = random.PRNGKey(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We keep our previousmy generated dataset but we will rely on `jnp.ndarray` instead of `np.ndarray`, by simply casting the arrays to JAX-compliant ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "x_jax = jnp.asarray(x_inputs)\n",
    "y_jax = jnp.asarray(y_class)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also keep the same hyperparameters as previously (see `eta`, `lambda_r` and `n_epochs` above) and randomly initialize our parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "eta = 0.01\n",
    "# Perform key splitting\n",
    "key_w, key_b, key = random.split(key, 3)\n",
    "# Initialize the weights\n",
    "weights = random.normal(key_w, (1, 2));\n",
    "bias = random.normal(key_b, (1, 1));\n",
    "# Save the weight history for plotting\n",
    "weights_history = jnp.zeros((n_epochs + 1, 2));\n",
    "bias_history = jnp.zeros((n_epochs + 1, 1));\n",
    "# Record the first \n",
    "weights_history = weights_history.at[0].set(weights[0])\n",
    "bias_history = bias_history.at[0].set(bias[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we discussed earlier, the original model by McCulloch & Pitts has several issues regarding its use for numerical optimization (as it was originally observed from a neuroscience point of view)\n",
    "\n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #410819; border-color: #cb2e47\">\n",
    "\n",
    "> **Changes from the `NumPy` implementation**\n",
    ">   1. (Smoothness) Our _activation function_ will change from _identity_ to sigmoid, i.e. $ \\phi(\\mathbf{x}) = \\frac{1}{1 + e^{-\\mathbf{x}}} $\n",
    ">   2. (Outliers) Our _loss function_ will change from $ L_{1} $ to $ L_{2} $ (MSE - Mean Squared Error)\n",
    "\n",
    "</div>\n",
    "\n",
    "***\n",
    "\n",
    "As discussed in the previous course, for using the very handy `grad` function provided by `JAX`, we need to define the whole forward function that we want to derive given our parameters. Hence this function also needs to **contain the forward pass**. This is summarized in the following function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import jax\n",
    "# Single neuron full function\n",
    "def loss_function(x, y, w, b):\n",
    "    # Forward pass\n",
    "    y_bar = jax.nn.sigmoid(jnp.dot(w, x.T) + b)\n",
    "    # Error computation\n",
    "    error = jnp.sum((y_bar - y) ** 2)\n",
    "    return error\n",
    "# Gradient of the loss\n",
    "grad_loss_function = jit(grad(loss_function, argnums=[2, 3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can now define the main training loop, which is almost exactly similar to the previously defined one minus the gradient computation operations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Update loop\n",
    "for i in range(n_epochs):\n",
    "\n",
    "    ######################\n",
    "    # Solution \n",
    "    gradients = grad_loss_function(x_jax, y_jax, weights, bias)\n",
    "    weights = weights - eta * gradients[0]\n",
    "    bias = bias - eta * gradients[1]\n",
    "    ######################\n",
    "    weights_history = weights_history.at[i + 1].set(weights[0])\n",
    "    bias_history = bias_history.at[i + 1].set(bias[0])\n",
    "    \n",
    "    if (verbose):\n",
    "        print('%2d. error = %f, weights = %f, %f, %f' % (i, np.sum(np.abs(error)) / n_observations, bias[0, 0], weights[0, 0], weights[0, 1]))\n",
    "plot = center_plot(\n",
    "    scatter_boundary(\n",
    "        np.array(x_inputs), \n",
    "        np.array(y_class), \n",
    "        np.array(weights_history), \n",
    "        np.array(bias_history), \n",
    "        i + 1)\n",
    "    )\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Multi-layer networks\n",
    "\n",
    "In the following, we define the overall architecture of the exercise to fill. The goal is to define a more advanced problem \n",
    "\n",
    "### 2-layer XOR problem\n",
    "\n",
    "In most cases, classification problems are far from being linear. Therefore, we need more advanced methods to be able to compute non-linear class boundaries. The advantage of neural networks is that the same principle can be applied in a *layer-wise* fashion. This allows to further discriminate the space in sub-regions (as seen in the course). We will try to implement the 2-layer *perceptron* that can provide a solution to the infamous XOR problem. The idea is now to have the output of the first neurons to be connected to a set of other neurons. Therefore, if we take back our previous formulation, we have the same output for the first neuron(s) $y$, that we will now term as $y^{(1)}$. Then, we feed these outputs to a second layer of neurons, which gives\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "y^{(2)}=\\sigma\\left(\\sum_{i = 1}^{n}w_{i}.y^{(1)}_{i} + b\\right)\n",
    "\\end{equation}\n",
    "$$  \n",
    "\n",
    "Finally, we will rely on the same loss $\\mathcal{L_{D}}$ as in the previous exercise, but the outputs used are $y^{(2)}$ instead of $y$. As in the previous case, we now need to compute the derivatives of the weights and biases for several layers . However, you should see that some form of generalization might be possible for any number of layer.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #192841; border-color: #779ecb\">\n",
    "\n",
    "> ### Exercise (course)\n",
    ">   1. Perform the derivatives for the last layer specifically\n",
    ">   2. Define a generalized derivative for any previous layer\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can construct the prototypical set of XOR values by using the following code (note that this is the most simple case, but still this is typically a problem that cannot be solved by a _linear classifier_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Input patterns\n",
    "x_inputs = np.array([[-1, -1],[-1,  1],[1, -1],[1,  1]])\n",
    "# Corresponding classes\n",
    "y_class = np.array([0, 1, 1, 0])                        \n",
    "# Initialize based on their sizes\n",
    "n_inputs = x_inputs.shape[0]\n",
    "n_outputs = 1\n",
    "# Plot the XOR problem\n",
    "plot = center_plot(scatter_classes(x_inputs, y_class, title = \"XOR classification\"))\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Observing the problem interactively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from cml.tasks import ClassificationLinear, ClassificationXOR\n",
    "explorer = ClassificationXOR()\n",
    "explorer.render()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "## Solving the problem with the `scikit-learn` library\n",
    "\n",
    "As a side note, and as discussed in the previous course, there is obviously a lot of libraries that can perform the training of such networks for you directly, without understanding what is exactly going on behind the scene. A quite extensive library is `scikit-learn`, which contains already coded models and learning procedure, that will allow us to _learn_ the parameters of this unknown function.\n",
    "\n",
    "Here we will use a `MLPClassifier` model to perfom classification on standardized features that we rescale through the `LinearRegression` and that this polynomial should be of degree 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPClassifier\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "# Create the MLP classifier\n",
    "classifier = MLPClassifier(alpha=1, max_iter=1000)\n",
    "# Standardize the input features\n",
    "clf = make_pipeline(StandardScaler(), classifier)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now based on this model definition, training it can be performed with a single line of code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from cml.data import xor_separation\n",
    "# Generate a XOR dataset\n",
    "x_inputs, y_classes = xor_separation(500, 0.1)\n",
    "# Train the MLP model\n",
    "clf.fit(x_inputs, y_classes);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Multi-layer implementation\n",
    "\n",
    "In this exercise, we will implement the XOR classification resolution using NumPy and JAX. This classification problem cannot be solved with a simple linear classifier. Hence, we will need to use at least a 2-layer network. The minimal set of variables that will be used by your code should be the following.\n",
    "\n",
    "```Python\n",
    "x_inputs          # 2 x n matrix of random points\n",
    "y_class           # classes of the patterns \n",
    "n_hidden          # Number of hidden units\n",
    "eta               # Learning rate parameter\n",
    "momentum          # Momentum parameter (bonus)\n",
    "weights1          # 1st layer weights\n",
    "weights2          # 2nd layer weights\n",
    "mse_limit         # Sum-squared error limit\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Number of hidden units\n",
    "n_hidden = 2\n",
    "# Learning rate parameter\n",
    "eta = 0.005\n",
    "# Momentum parameter\n",
    "momentum = 0.1\n",
    "# Sum-squared error limit\n",
    "mse_limit = 0.002"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Exercise 1 - 2-layers XOR classification\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #013220; border-color: #03C03C\">\n",
    "\n",
    "> To complete this exercise, you will need to have a basic understanding of NumPy and JAX. You can use the resources provided in the previous exercises to learn more about these libraries.\n",
    "\n",
    "We help you out by first defining the XOR classification problem to solve (and outlining some basic properties)\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from cml.data import xor_separation\n",
    "# Definition of our data properties\n",
    "n_observations = 500\n",
    "noise_factor = 0.1\n",
    "# Properties of our problem\n",
    "n_input = 2\n",
    "n_output = 1\n",
    "# Generate a XOR dataset\n",
    "x_inputs, y_classes = xor_separation(n_observations, noise_factor)\n",
    "\n",
    "for i in range(y_classes.shape[0]):\n",
    "    if(y_classes[i] == 0):\n",
    "        y_classes = y_classes.at[i].set(-1)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also show some potential implementation of the **initialization scheme**. Note that you can later try to replace this naive implementation by some more refined functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 1st layer weights\n",
    "weights1 = (np.random.randn(n_hidden, n_input) * 0.01)\n",
    "# 2nd layer weights\n",
    "weights2 = (np.random.randn(n_output, n_hidden) * 0.01)\n",
    "# 1st layer biases\n",
    "bias1 = (np.random.randn(n_hidden, 1) * 0.01)\n",
    "# 2nd layer biases\n",
    "bias2 = (np.random.randn(n_output, 1) * 0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #013220; border-color: #03C03C\">\n",
    "\n",
    "> ### Question 1.1 - Implement 2-layer XOR classification with NumPy\n",
    "\n",
    "> 1. Update the forward propagation and error computation\n",
    "> 2. Update the back-propagation part to learn the weights of both layers.\n",
    "> 3. Run the learning and check your network results\n",
    "  \n",
    "*For optional questions, please look at the end of this exercise's code boxes for more information*\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#1. Calculating forward propagation\n",
    "n_iter = 10\n",
    "\n",
    "\n",
    "def sigmoid(x):\n",
    "    return 1/(1 + np.exp(-x))\n",
    "\n",
    "def hinge_loss(y, y_bar):\n",
    "    return max(0,1-y*y_bar)\n",
    "\n",
    "\n",
    "for _ in range(n_iter):\n",
    "    #forward propagation\n",
    "    o1 = np.matmul(weights1,x_inputs.transpose())\n",
    "    y1 = sigmoid(o1)#this is the output after passing through the first hidden node\n",
    "\n",
    "    o2 = np.matmul(weights2,y1)\n",
    "    y2 = sigmoid(o2)\n",
    "\n",
    "    y2 = y2.transpose()\n",
    "    y1 = y1.transpose()\n",
    "    o1 = o1.transpose()\n",
    "    o2 = o2.transpose()\n",
    "\n",
    "    #error computation using Hinge Loss\n",
    "    #Back-propagation to learn the weights of both layers\n",
    "    dw2_11 = 0\n",
    "    dw2_12 = 0\n",
    "    dw1_11 = 0\n",
    "    dw1_12 = 0\n",
    "    dw1_21 = 0\n",
    "    dw1_22 = 0\n",
    "    \n",
    "    error = 0\n",
    "    for i in range(y2.shape[0]):\n",
    "        error += hinge_loss(y_classes[i],y2[i])\n",
    "\n",
    "    for i in range(y2.shape[0]):\n",
    "\n",
    "        if(y_classes[i]*y2[i] < 1):\n",
    "            #calculation of necessary derivatives for backpropagation\n",
    "            \n",
    "            #second layer of weights\n",
    "            dL_dy2 = -y_classes[i]*y2[i]\n",
    "            dy2_do2 = sigmoid(o2[i])/(1-sigmoid(o2[i]))\n",
    "\n",
    "            do2_dw2_11 = y1[i][0]\n",
    "            do2_dw2_21 = y1[i][1]\n",
    "\n",
    "            dL_dw2_11 = dL_dy2*dy2_do2*do2_dw2_11\n",
    "            dL_dw2_21 = dL_dy2*dy2_do2*do2_dw2_21\n",
    "\n",
    "            #first layer of weights\n",
    "            do2_dy1_1 = weights2[0][0]\n",
    "            dy1_1_do_11 = sigmoid(o1[i][0])/(1-sigmoid(o1[i][0]))\n",
    "            do_11_dw1_11 = x_inputs[i][0] \n",
    "            dL_dw1_11 = dL_dy2*dy2_do2*do2_dy1_1*dy1_1_do_11*do_11_dw1_11\n",
    "\n",
    "            do_11_dw1_12 = x_inputs[i][1]\n",
    "            dL_dw1_12 = dL_dy2*dy2_do2*do2_dy1_1*dy1_1_do_11*do_11_dw1_12\n",
    "\n",
    "            do2_dy2_1 = weights2[0][1]\n",
    "            dy2_1_do_12 = sigmoid(o1[i][1])/(1-sigmoid(o1[i][1]))\n",
    "            do_12_dw1_21 = x_inputs[i][0] \n",
    "            dL_dw1_21  = dL_dy2*dy2_do2*do2_dy2_1*dy2_1_do_12*do_12_dw1_21\n",
    "\n",
    "            do_12_dw1_22 = x_inputs[i][1] \n",
    "            dL_dw1_22  = dL_dy2*dy2_do2*do2_dy2_1*dy2_1_do_12*do_12_dw1_22\n",
    "\n",
    "            dw2_11 += dL_dw2_11\n",
    "            dw2_12 += dL_dw2_12\n",
    "            dw1_11 += dL_dw1_11\n",
    "            dw1_21 += dL_dw1_21\n",
    "            dw1_12 += dL_dw1_12\n",
    "            dw1_22 += dL_dw1_22\n",
    "\n",
    "    weights1[0][0] -= eta*dw1_11\n",
    "    weights1[0][1] -= eta*dw1_12\n",
    "    weights1[1][0] -= eta*dw1_21\n",
    "    weights1[1][1] -= eta*dw1_22\n",
    "    weights2[0][0] -= eta*dw2_11\n",
    "    weights2[0][1] -= eta*dw2_12\n",
    "\n",
    "    print(weights1)\n",
    "    print(weights2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #013220; border-color: #03C03C\">\n",
    "\n",
    "> ### Question 1.2 - Implement XOR classification with JAX\n",
    "\n",
    "> 1. Perform the same implementation as previously with JAX\n",
    "> 2. Implement gradient descent to optimize the weights and bias.\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_iter = 10\n",
    "\n",
    "\n",
    "def sigmoid(x):\n",
    "    return 1/(1 + jnp.exp(-x))\n",
    "\n",
    "def hinge_loss_jnp(y, x_inputs, weights1, weights2):\n",
    "    o1 = jnp.matmul(weights1,x_inputs.transpose())\n",
    "    y1 = sigmoid(o1)\n",
    "    o2 = jnp.matmul(weights2,y1)\n",
    "    y2 = sigmoid(o2)\n",
    "    return jnp.maximum(0,1-(y*y2)).sum() #sum() just converts datatype. it's already of size 1\n",
    "                    \n",
    "grad_loss_function = jit(grad(hinge_loss_jnp, argnums=[2, 3]))\n",
    "\n",
    "for _ in range(n_iter):\n",
    "    #forward propagation\n",
    "   \n",
    "    dw2_11 = 0\n",
    "    dw2_12 = 0\n",
    "    dw1_11 = 0\n",
    "    dw1_12 = 0\n",
    "    dw1_21 = 0\n",
    "    dw1_22 = 0\n",
    "    \n",
    "    for i in range(x_inputs.shape[0]):\n",
    "\n",
    "        #calculation of necessary derivatives for backpropagation\n",
    "\n",
    "        gradients = grad_loss_function(y_classes[i],x_inputs[i],weights1,weights2)\n",
    "                \n",
    "        dw2_11 += gradients[1][0][0]\n",
    "        dw2_12 += gradients[1][0][1]\n",
    "        dw1_11 += gradients[0][0][0]\n",
    "        dw1_21 += gradients[0][1][0]\n",
    "        dw1_12 += gradients[0][0][1]\n",
    "        dw1_22 += gradients[0][1][1]\n",
    "\n",
    "    weights1[0][0] -= eta*dw1_11\n",
    "    weights1[0][1] -= eta*dw1_12\n",
    "    weights1[1][0] -= eta*dw1_21\n",
    "    weights1[1][1] -= eta*dw1_22\n",
    "    weights2[0][0] -= eta*dw2_11\n",
    "    weights2[0][1] -= eta*dw2_12\n",
    "\n",
    "    print(weights1)\n",
    "    print(weights2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #013220; border-color: #03C03C\">\n",
    "\n",
    "> ### Question 1.3 - Understanding properties of our implementation\n",
    "\n",
    "> 3. Perform multiple re-runs of the learning procedure (re-launching with different initializations)\n",
    ">     1. What observations can you make on the learning process?\n",
    ">     2. What happens if you initialize all weights to zeros?\n",
    ">     3. Change your initialization and regularization scheme\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "## Question 1.3: Answers\n",
    "# 3. Perform multiple re-runs of the learning procedure (re-launching with different initializations)\n",
    "#     1. Depending on the chosen initializations (iteration number, learning rate, weight initialization, the achieved accuracy varies).\n",
    "#     2. Initalizing all weights with 0 will mean that no values passed forward on the first iteration. The outputs will therefore all be 0, which will not pose an error, so gradient descent will still lead to a change weights.\n",
    "#     3. Code below for changing the regularization and initialization. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Using the jax implementation:\n",
    "\n",
    "n_iter = 10\n",
    "\n",
    "#regularization parameters\n",
    "l2 = True\n",
    "l1 = True\n",
    "\n",
    "#regularization factor\n",
    "lamb = 0.01 \n",
    "\n",
    "def sigmoid(x):\n",
    "    return 1/(1 + jnp.exp(-x))\n",
    "\n",
    "#loss function\n",
    "def hinge_loss_jnp(y, x_inputs, weights1, weights2):\n",
    "    o1 = jnp.matmul(weights1,x_inputs.transpose())\n",
    "    y1 = sigmoid(o1)\n",
    "    o2 = jnp.matmul(weights2,y1)\n",
    "    y2 = sigmoid(o2)\n",
    "    return jnp.maximum(0,1-(y*y2)).sum() #sum() just converts datatype. it's already of size 1\n",
    "                 \n",
    "grad_loss_function = jit(grad(hinge_loss_jnp, argnums=[2, 3]))\n",
    "\n",
    "#regularization selection \n",
    "if(l1 and not l2):\n",
    "    lambl1 = lamb\n",
    "    lambl2 = 0\n",
    "elif(l2 and not l1):\n",
    "    lambl1 = 0\n",
    "    lambl2 = lamb\n",
    "elif(l1 and l2):\n",
    "    lambl1 = lamb\n",
    "    lambl2 = lamb\n",
    "else:\n",
    "    lambl1 = 0\n",
    "    lambl2 = 0\n",
    "\n",
    "for _ in range(n_iter):\n",
    "    \n",
    "    dw2_11 = 0\n",
    "    dw2_12 = 0\n",
    "    dw1_11 = 0\n",
    "    dw1_12 = 0\n",
    "    dw1_21 = 0\n",
    "    dw1_22 = 0\n",
    "    \n",
    "    for i in range(x_inputs.shape[0]):\n",
    "\n",
    "        #calculation of necessary derivatives for backpropagation\n",
    "\n",
    "        gradients = grad_loss_function(y_classes[i],x_inputs[i],weights1,weights2)\n",
    "               \n",
    "        dw2_11 += gradients[1][0][0]\n",
    "        dw2_12 += gradients[1][0][1]\n",
    "        dw1_11 += gradients[0][0][0]\n",
    "        dw1_21 += gradients[0][1][0]\n",
    "        dw1_12 += gradients[0][0][1]\n",
    "        dw1_22 += gradients[0][1][1]\n",
    "    \n",
    "    #adding the regularization\n",
    "\n",
    "    weights1[0][0] -= eta*(dw1_11 + lambl1 + lambl2*weights1[0][0])\n",
    "    weights1[0][1] -= eta*(dw1_12 + lambl1 + lambl2*weights1[0][1])\n",
    "    weights1[1][0] -= eta*(dw1_21 + lambl1 + lambl2*weights1[1][0])\n",
    "    weights1[1][1] -= eta*(dw1_22 + lambl1 + lambl2*weights1[1][1])\n",
    "    weights2[0][0] -= eta*(dw2_11 + lambl1 + lambl2*weights2[0][0])\n",
    "    weights2[0][1] -= eta*(dw2_12 + lambl1 + lambl2*weights2[0][1])\n",
    "\n",
    "    print(weights1)\n",
    "    print(weights2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #192841; border-color: #779ecb\">\n",
    "\n",
    "> ### Going further\n",
    "\n",
    "> The following set of questions are optionnal addendum to your previous code that allow to understand more in-depth properties about _regularization_ both for neural networks and optimization in general.\n",
    "\n",
    "> 1. (Optional) Implement the *weight decay* constraint in your network.\n",
    "> 2. (Optional) Add the *momentum* to the learning procedure.\n",
    "\n",
    "> *Weight decay* constraint\n",
    "> As nothing constrains the weights in the network, we can note that usually all weights vector given a multiplicative factor might be equivalent, which can stall the learning (and lead to exploding weights). The *weight decay* allows to regularize the learning by penalizing weights with a too wide amplitude. The idea is to add this constraint as a term to the final loss (which leads to an indirect \"pressure\" on the learning process. Therefore, the final loss will be defined as\n",
    "> $$\n",
    "\\begin{equation}\n",
    "\\mathcal{L}_{final}=\\mathcal{L_D} + \\lambda \\sum_{l} \\sum_{i} \\sum_{j} \\left( W_{ij}^{l} \\right)^{2}\n",
    "\\end{equation}\n",
    "$$\n",
    "> where the parameter $\\lambda$ controls the relative importance of the two terms.\n",
    "\n",
    "> *Momentum* in learning\n",
    "> Usually, in complex problems, the gradient can be very noisy and, therefore, the learning might oscillate widely. In order to reduce this problem, we can *smooth* the different gradient updates by retaining the values of the gradient at each iteration and then performing an update based on the latest gradient $\\delta_{i}^{t}$ and the gradient at the previous iteration $\\delta_{i}^{t-1}$. Therefore, a gradient update is applied as\n",
    "> $$\n",
    "\\begin{equation}\n",
    "\\delta_{final}^{t} = \\delta_{i}^{t} + m.\\delta_{i}^{t-1}\n",
    "\\end{equation}\n",
    "$$\n",
    "> with $m$ the momentum parameter, which control the amount of gradient smoothing.\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "######################\n",
    "# YOUR CODE GOES HERE\n",
    "######################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (Optional) 3-layer audio classification\n",
    "\n",
    "Once again, note that the following paragraph and subsequent questions are _optional_ but should be a quite simple extension of our previous work on 2-layer networks. If you struggle with the data import mechanisms, you can look down at the PyTorch section (with a **mandatory exercise**) that will provide base code for this aspect. \n",
    "\n",
    "Finally, we will attack a complete audio classification problem and try to perform neural network learning on a set of audio files. The data structure will be the same as the one used for parts 1 and 2. As discussed during the courses, even though a 2-layer neural network can provide non-linear boundaries, it can not perform \"holes\" inside those regions. In order to obtain an improved classification, we will now rely on a 3-layer neural network. The modification to the code of section 3.2 should be minimal, as the back-propagation will be similar for the new layer as one of the two others. We do not develop the math here as it is simply a re-application of the previous rules with an additional layer (which derivatives you should have generalized in the previous exercise).  \n",
    "\n",
    "However, up until now, we only performed *binary classification* problems, but this time we need to obtain a decision rule for multiple classes. Therefore, we cannot rely on simply computing the distance between desired patterns and the obtained binary value. The idea here is to rely on the *softmax regression*, by considering classes as a vector of probabilities. The desired answers will therefore be considered as a set of *probabilities*, where the desired class is $1$ and the others are $0$ (called *one-hot* representation). Then, the cost function will rely on the softmax formulation\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathcal{L_D}(\\theta) = - \\frac{1}{m} \\left[ \\sum_{i=1}^{m} \\sum_{j=1}^{k} 1\\left\\{y^{(i)} = j\\right\\} \\log \\frac{e^{\\theta_{j}^{T} x^{(i)}}}{\\sum_{l=1}^{k} e^{ \\theta_{l}^{T} x^{(i)} }}  \\right]\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Therefore, we compute the output of the softmax by taking \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "p(y^{(i)} = j | x^{(i)}; \\theta) = \\frac{e^{\\theta_{j}^{T} x^{(i)}}}{\\sum_{l=1}^{k} e^{ \\theta_{l}^{T} x^{(i)}} }\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "By taking derivatives, we can show that the gradient of the softmax layer is\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\nabla_{\\theta_{j}} \\mathcal{L_D}(\\theta) = - \\frac{1}{m} \\sum_{i=1}^{m}{ \\left[ x^{(i)} \\left( 1\\{ y^{(i)} = j\\}  - p(y^{(i)} = j \\mid x^{(i)}, \\theta) \\right) \\right]}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sweet activation functions\n",
    "\n",
    "As discussed in the course, the interest of stacking layers is that there is an _activation function_, which allows non-linear interactions between the dimensions (and avoids to only compute a single huge affine transform). Although the `sigmoid` function has been historically the most used, there has been some large developments since. Notably the `ReLU` (Rectified Linear Unit) is one of the major difference in modern networks (we will see more about that in a later course)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function for computing the Sigmoid activation\n",
    "def sigmoid(x):\n",
    "    return 1 / (1 + np.exp(-x))\n",
    "# Derivative\n",
    "def dsigmoid(a):\n",
    "    return a * (1.0 - a)\n",
    "# Function for computing the ReLU activation\n",
    "def relu(x):\n",
    "    return np.maximum(0, x)\n",
    "# Derivative\n",
    "def drelu(x):\n",
    "    if (x < 0):\n",
    "        return 0\n",
    "    return 1\n",
    "# Function for computing the Tanh activation\n",
    "def tanh(x):\n",
    "    return np.tanh(x);\n",
    "# Derivative\n",
    "def dtanh(x): \n",
    "    return np.cosh(x) ^ -2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we plot some simple examples of what these activation functions look like. You can try to rely on these functions in your previous training code and witness the differences in training.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.models import Div\n",
    "from bokeh.layouts import column, row\n",
    "from cml.plot import cml_figure\n",
    "# Functions to\n",
    "funcs = [('Sigmoid',sigmoid,'red'), ('Tanh',tanh,'orange'), ('ReLU',relu,'yellow')]\n",
    "# Generating the x axis\n",
    "x = np.linspace(-5, 5, 100)\n",
    "plots = []\n",
    "for (name, func, color) in funcs:\n",
    "    cur_plot = cml_figure(plot_width=400, plot_height=250, title=name)\n",
    "    cur_plot.line(x, func(x), color=color, line_width=4)\n",
    "    plots.append(cur_plot)\n",
    "plot = center_plot(column(Div(text = \"Different activation functions\", style={'font-size': '250%'}), row(*plots)))\n",
    "plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coding the whole network from scratch\n",
    "\n",
    "You should now have all the tools necessary to apply neural networks from scratch to a more complex problem. In the following exercise, we simply removed any guideline code, and you need to code all the procedure for training a NN and **apply it to audio data**. You will use the spectral features discussed in the previous exercise as an input.\n",
    "\n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #192841; border-color: #779ecb\">\n",
    "\n",
    "> ### Going further **(optional exercise) **\n",
    "\n",
    ">  1. Based on the previous neural network, upgrade the code to a 3-layer neural network\n",
    ">  2. Implement the *softmax regression* on top of your 3-layer network\n",
    ">  3. Use the provided code to perform classification on a pre-defined set of features\n",
    ">  4. As previously, change the set of features to assess their different accuracies\n",
    ">  5. Evaluate the neural network accuracy for all features combinations\n",
    ">  6. What happens if the learning rate is too large ? What is this phenomenon ?\n",
    ">  7. Perform a more advanced visualization of the learning process.\n",
    "  \n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "######################\n",
    "# YOUR CODE GOES HERE\n",
    "######################\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using `Pytorch` to enjoy life\n",
    "\n",
    "Up to now, we have been writing every operations by ourselves (in order to better understand the mathematics behind NN). However, there exists of course some simplifying libraries that provide large simplifications to this question.\n",
    "\n",
    "One of the most powerful and complete library of this sort is `Pytorch`, which has been developed for several years (even prior to the recent boom of deep learning). `Pytorch` provides a large set of pre-coded layers, but also **computational graphs** and **autograd**, which are very powerful paradigms allowing to define complex operators and automatically taking derivatives.\n",
    "\n",
    "## An (extremely) fast and dirty introduction to `Pytorch`\n",
    "\n",
    "PyTorch is a popular open-source deep learning framework based on the Torch library. It is primarily developed by Facebook AI research team, and it provides a simple and efficient way to build deep learning models. In this extremely short crash course, we will cover the basics of PyTorch, including:\n",
    "\n",
    "* Tensors\n",
    "* Automatic differentiation\n",
    "* Neural networks\n",
    "* Training a neural network on a dataset\n",
    "\n",
    "For a complete tutorial on PyTorch, please check out the [official PyTorch documentation](https://pytorch.org/tutorials/). This tutorial covers a wide range of topics, including the basics of PyTorch, building and training neural networks, working with data and dataloaders, and deploying models to production."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tensors\n",
    "At its core, PyTorch is all about tensors. A tensor is a generalization of vectors and matrices to an arbitrary number of dimensions. A scalar is a 0-dimensional tensor, a vector is a 1-dimensional tensor, and a matrix is a 2-dimensional tensor.\n",
    "\n",
    "#### Creating Tensors\n",
    "We can create a PyTorch tensor from a Python list or a NumPy array using the `torch.tensor()` function, but also tensors of a specific size and/or with all zeros or all ones using the `torch.zeros()` and `torch.ones()` functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "# create a tensor from a Python list\n",
    "x = torch.tensor([1, 2, 3, 4])\n",
    "print(x)\n",
    "# create a tensor from a NumPy array\n",
    "y = torch.tensor(np.array([[1, 2], [3, 4]]))\n",
    "print(y)\n",
    "# create a 2x3 tensor of zeros\n",
    "z = torch.zeros((2, 3))\n",
    "print(z)\n",
    "# create a 3x2 tensor of ones\n",
    "w = torch.ones((3, 2))\n",
    "print(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tensor operations\n",
    "\n",
    "PyTorch tensors support a wide range of operations, including arithmetic operations like addition, subtraction, multiplication, and division, as well as more advanced operations like matrix multiplication, element-wise multiplication, and concatenation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create two tensors\n",
    "x = torch.tensor([1, 2, 3, 4])\n",
    "y = torch.tensor([5, 6, 7, 8])\n",
    "# add the two tensors\n",
    "z = x + y\n",
    "print(z)\n",
    "# element-wise multiplication\n",
    "w = x * y\n",
    "print(w)\n",
    "# matrix multiplication\n",
    "a = torch.tensor([[1, 2], [3, 4]])\n",
    "b = torch.tensor([[5, 6], [7, 8]])\n",
    "c = torch.matmul(a, b)\n",
    "print(c)\n",
    "# concatenation\n",
    "d = torch.cat((a, b), dim=1)\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic differentiation\n",
    "\n",
    "The core of PyTorch is the autograd package. It provides automatic differentiation for all operations on Tensors.\n",
    "\n",
    "PyTorch's autograd package provides automatic differentiation for all operations on Tensors. It is a define-by-run framework, which means that your backpropagation is defined by how your code is run, and that every single iteration can be different.\n",
    "\n",
    "To compute gradients, the `requires_grad` property of Tensors needs to `True` (which is the case by default). This tells PyTorch to track all operations on the Tensor. You can then call `.backward()` on a Tensor to compute the gradients with respect to that Tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "x = torch.tensor(2.0, requires_grad=True)\n",
    "y = 2*x**2 + 3*x - 1\n",
    "y.backward()\n",
    "\n",
    "print(x.grad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining our network\n",
    "\n",
    "When building neural networks we frequently think of arranging the computation into layers, some of which have learnable parameters which will be optimized during learning. In `PyTorch`, the `nn` package provides higher-level abstractions over raw computational graphs that are useful for building neural networks. The `nn` package defines a set of `Modules`, which are roughly equivalent to neural network layers. A `Module` receives input `Tensors` and computes output `Tensors`, but may also hold internal state such as `Tensors` containing learnable parameters. The nn package also defines a set of useful loss functions that are commonly used when training neural networks.\n",
    "\n",
    "In the following example, we use the `nn` package to show how easy it is to instantiate our previous three-layers network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "# Define the input dimensions\n",
    "in_size = 1000\n",
    "# Number of neurons in a layer\n",
    "hidden_size = 100\n",
    "# Output (target) dimension\n",
    "output_size = 10\n",
    "# Use the nn package to define our model and loss function.\n",
    "model = torch.nn.Sequential(\n",
    "    torch.nn.Linear(in_size, hidden_size),\n",
    "    torch.nn.Sigmoid(),\n",
    "    torch.nn.Linear(hidden_size, hidden_size),\n",
    "    torch.nn.Tanh(),\n",
    "    torch.nn.Linear(hidden_size, output_size),\n",
    "    torch.nn.Softmax(dim=1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimizing the network\n",
    "\n",
    "Up to this point we have updated the weights of our models by manually performing the gradient descent algorithm (changing the parameters vectors). Although this is not a huge burden for simple optimization algorithms like stochastic gradient descent, in practice we often train neural networks using more sophisticated optimizers like AdaGrad, RMSProp or Adam (that we will see later in this course)\n",
    "\n",
    "The `optim` package in PyTorch abstracts the idea of an optimization algorithm and provides implementations of commonly used optimization algorithms, and greatly simplfies the training loop associated with training a neural network.\n",
    "\n",
    "For the sake of presentation we will use random inputs $\\mathbf{x}$ that should be matched with random outputs $\\mathbf{y}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_size = 64\n",
    "# Create random Tensors to hold inputs and outputs\n",
    "x = torch.randn(batch_size, in_size)\n",
    "y = torch.randn(batch_size, output_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following example we optimize the model using the Adam algorithm provided by the `optim` package, based on a `MSE` loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learning rate\n",
    "learning_rate = 1e-4\n",
    "# Loss function that we will use\n",
    "loss_fn = torch.nn.MSELoss(reduction='sum')\n",
    "# Optimizer to fit the weights of the network\n",
    "optimizer = torch.optim.SGD(model.parameters(), lr=learning_rate)\n",
    "for t in range(500):\n",
    "    # Forward pass: compute predicted y by passing x to the model.\n",
    "    y_pred = model(x)\n",
    "    # Compute the loss.\n",
    "    loss = loss_fn(y_pred, y)\n",
    "    # Before the backward pass, zero all of the network gradients\n",
    "    optimizer.zero_grad()\n",
    "    # Backward pass: compute gradient of the loss with respect to parameters\n",
    "    loss.backward()\n",
    "    # Calling the step function to update the parameters\n",
    "    optimizer.step()\n",
    "print(loss)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 - Pytorch to classify audio\n",
    "\n",
    "Now that we know the main components of `Pytorch` to define and optimize networks, your assignement is to define a complete classification problem from audio data, by relying on this toolbox\n",
    "\n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #013220; border-color: #03C03C\">\n",
    "\n",
    "> ### Question 1.1 - Model and definition\n",
    "\n",
    ">   1. Use `Pytorch` to define a model for audio classification\n",
    ">   2. Import the audio features dataset and check that your model produces an output\n",
    "\n",
    "</div>\n",
    "\n",
    "***\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cml.data import import_dataset, AudioSupervisedDataset\n",
    "data_path = \"/home/frederik/creative_ml/data\"\n",
    "batch_size = 1\n",
    "# Instantiate the dataset class\n",
    "train_dataset = import_dataset(data_path, \"musclefish\", 'train', batch_size=batch_size)\n",
    "# Generate a batch of data from the train dataset\n",
    "batch_x, batch_y = next(train_dataset)\n",
    "# Print the shape of the input and output batch tensors\n",
    "print(batch_x.shape)  # (16, 44100)\n",
    "print(batch_y.shape)  # (16,)\n",
    "#dataset has 44100 different sampling points per input and we segment them in batches of 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# power_spec = train_dataset.transform(10, \"stft\")\n",
    "# print(power_spec.shape)\n",
    "\n",
    "import torch\n",
    "import numpy as np\n",
    "\n",
    "# Define the input dimensions\n",
    "in_size = 44100\n",
    "# Number of neurons in a layer\n",
    "hidden_size = 20\n",
    "# Output (target) dimension\n",
    "output_size = 16\n",
    "# Use the nn package to define our model and loss function.\n",
    "\n",
    "model = torch.nn.Sequential(\n",
    "    torch.nn.Linear(in_size, hidden_size),\n",
    "    torch.nn.Sigmoid(),\n",
    "    torch.nn.Linear(hidden_size, hidden_size),\n",
    "    torch.nn.Sigmoid(),\n",
    "    torch.nn.Linear(hidden_size, hidden_size),\n",
    "    torch.nn.Sigmoid(),\n",
    "    torch.nn.Linear(hidden_size, output_size),\n",
    "    torch.nn.Softmax(dim=1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "***\n",
    "\n",
    "<div class=\"alert alert-success\" markdown=\"1\" style=\"color:white; background-color: #013220; border-color: #03C03C\">\n",
    "\n",
    "> ### Question 1.2 - Training the model\n",
    "\n",
    ">   3. Write the optimization loop (think carefully about the _loss function_\n",
    ">   4. As previously, change the set of features to assess their different accuracies\n",
    ">   5. (Optional) Think of how you could use more complex features (time series, audio, STFT) to classify your data\n",
    "\n",
    "</div>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Learning rate\n",
    "learning_rate = 1e-4\n",
    "\n",
    "# Loss function that we will use\n",
    "loss_fn = torch.nn.CrossEntropyLoss()\n",
    "\n",
    "# Optimizer to fit the weights of the network\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate)\n",
    "\n",
    "train_dataset._reset_generator()\n",
    "for batch_x, batch_y in iter(train_dataset):\n",
    "    \n",
    "    torch_x = torch.from_numpy(batch_x.numpy())\n",
    "    torch_y = torch.from_numpy(batch_y.numpy())\n",
    "    torch_y = torch_y.long()\n",
    "    \n",
    "    y_pred = model(torch_x)\n",
    "    loss = loss_fn(y_pred, torch_y)\n",
    "    optimizer.zero_grad()\n",
    "    loss.backward()\n",
    "\n",
    "    optimizer.step()\n",
    "    \n",
    "print(loss)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#4.\n",
    "lr_vector = np.linspace(1e-6,5e-2,30)\n",
    "loss_lr = np.zeros((30,1))\n",
    "\n",
    "for i,lr in enumerate(lr_vector):\n",
    "    train_dataset._reset_generator()\n",
    "    for batch_x, batch_y in iter(train_dataset):\n",
    "    \n",
    "        torch_x = torch.from_numpy(batch_x.numpy())\n",
    "        torch_y = torch.from_numpy(batch_y.numpy())\n",
    "        torch_y = torch_y.long()\n",
    "\n",
    "        y_pred = model(torch_x)\n",
    "        loss = loss_fn(y_pred, torch_y)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "\n",
    "        optimizer.step()\n",
    "    loss_lr[i] = loss.detach().numpy()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(lr_vector,loss_lr,'ro')\n",
    "plt.xlabel('Learning rate')\n",
    "plt.ylabel('Absolute loss')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print('We can observe that the learning rate leads to a significant jump in loss; however, there does not seem to be one single minimum, but rather an accumulation of local minima.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}