{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Multilayer Perceptron\n\n\nIn this notebook, we'll start with one of the most basic neural network architectures, a multilayer perceptron (MLP), also known as a feedforward network. The dataset we'll be using is the famous MNIST dataset, a dataset of 28x28 black and white images consisting of handwritten digits, 0 to 9.\n\n![](https://github.com/bentrevett/pytorch-image-classification/blob/master/assets/mlp-mnist.png?raw=1)\n\nWe'll process the dataset, build our model and then train our model. Afterwards we'll do a short dive into what the model has actually learned.\n\n### Data Processing\n\nLet's start by importing all of the modules we'll need. The main ones we need to import are:\n- torch for general PyTorch functionality\n- torch.nn and torch.nn.functional for neural network based functions\n- torch.optim for our optimizer which will update the parameters of our neural network\n- torch.utils.data for handling the dataset\n- torchvision.transforms for data augmentation\n- torchvision.datasets for loading the dataset\n- sklearn's metrics for visualizing a confusion matrix\n- sklearn's decomposition and manifold for visualizing the neural network's representations in two dimensions\n- matplotlib for plotting"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim\nimport torch.utils.data as data\n\nimport torchvision.transforms as transforms\nimport torchvision.datasets as datasets\n\nfrom sklearn import metrics\nfrom sklearn import decomposition\nfrom sklearn import manifold\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport copy\nimport random\nimport time","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To ensure we get reproducible results we set the random seed for Python, Numpy and PyTorch."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"SEED = 1234\n\nrandom.seed(SEED)\nnp.random.seed(SEED)\ntorch.manual_seed(SEED)\ntorch.cuda.manual_seed(SEED)\ntorch.backends.cudnn.deterministic = True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first thing we'll do is load our dataset.\n\nThis will automatically download the training set for the MNIST dataset and save it in a folder called `.data`. It will create the folder if it does not exist."},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"ROOT = '.data'\n\ntrain_data = datasets.MNIST(root = ROOT, \n                            train = True, \n                            download = True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we want to *normalize* our data. This means we want it to have a mean of zero and a standard deviation of one. \n\nWhy do we want to do this? Normalizing our data allows our models to train faster and to also help them avoid local minima, i.e. train more reliably.\n\nWe normalize our data by subtracting the mean and dividing by the standard deviation of our dataset. First, we need to calculate the mean and standard deviation. <br>\n**Note**: It is important that the mean and standard deviation are only calculated over the training set and not the test set. We do not want to use any information at all from the test set and should only look at it when we are calculating our test loss.\n\nTo calculate the means and standard deviations we get the actual data (the images) using the `.data.` attribute of our training data, convert them into floating point numbers and then use the built in `mean` and `std` functions to calculate the mean and standard deviation, respectively. The image data has values between 0-255, which we want to scale between 0-1, so we divide by 255."},{"metadata":{"trusted":true},"cell_type":"code","source":"mean = train_data.data.float().mean() / 255\nstd = train_data.data.float().std() / 255","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Calculated mean: {mean}')\nprint(f'Calculated std: {std}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing The Data\nNow we've calculated our mean and standard deviation how do we actually use them? We use Torchvision's `transforms`. \n\nA `transform` states how our data should be augmented and processed. Data augmentation involves manipulating the available training data in a way that artifically creates more training examples. We use `transforms.Compose` to built a list of transformations that will be applied to the image. \n\nThe transforms we use are:\n- `RandomRotation` - randomly rotates the image between `(-x, +x)` degrees, where we have set `x = 5`. Note the `fill=(0,)` is due to a [bug](https://github.com/pytorch/vision/issues/1759) in some versions of torchvision. \n- `RandomCrop` - this first adds `padding` around our image, 2 pixels here, to artifically make it bigger, before taking a random `28x28` square crop of the image.\n- `ToTensor()` - this converts the image from a PIL image into a PyTorch tensor.\n- `Normalize` - this subtracts the mean and divides by the standard deviations given. \n\nThe first two transformations have to be applied before `ToTensor` as they should both be applied on a PIL image. `Normalize` should only be applied to the images after they have been converted into a tensor. See the Torchvision documentation for [transforms that should be applied to PIL images](https://pytorch.org/docs/stable/torchvision/transforms.html#transforms-on-pil-image) and [transforms that should be applied on tensors](https://pytorch.org/docs/stable/torchvision/transforms.html#transforms-on-torch-tensor).\n\nWe have two lists of transforms, a train and a test transform. The train transforms are to artificially create more examples for our model to train on. We do not augment our test data in the same way as we want a consistent set of examples to evaluate our final model on. The test data, however, should still be normalized."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_transforms = transforms.Compose([\n                            transforms.RandomRotation(5, fill=(0,)),\n                            transforms.RandomCrop(28, padding = 2),\n                            transforms.ToTensor(),\n                            transforms.Normalize(mean = [mean], std = [std])\n                                      ])\n\ntest_transforms = transforms.Compose([\n                           transforms.ToTensor(),\n                           transforms.Normalize(mean = [mean], std = [std])\n                                     ])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now we have defined our transforms we can then load the train and test data with the relevant transforms defined."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data = datasets.MNIST(root = ROOT, \n                            train = True, \n                            download = True, \n                            transform = train_transforms)\n\ntest_data = datasets.MNIST(root = ROOT, \n                           train = False, \n                           download = True, \n                           transform = test_transforms)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can simply check the `len` of the datasets to see how many examples are within each."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Number of training examples: {len(train_data)}')\nprint(f'Number of testing examples: {len(test_data)}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can get a look at some of the images within our dataset to see what we're working with. The function below plots a square grid of images. If you supply less than a complete square number of images it will ignore the last few."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_images(images):\n\n    n_images = len(images)\n\n    rows = int(np.sqrt(n_images))\n    cols = int(np.sqrt(n_images))\n\n    fig = plt.figure()\n    for i in range(rows*cols):\n        ax = fig.add_subplot(rows, cols, i+1)\n        ax.imshow(images[i].view(28, 28).cpu().numpy(), cmap = 'bone')\n        ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's load 25 images. These will have been processed through our transforms so will be randomly rotated and cropped.\n\nIt's good practice to see your data with your transforms applied so you can ensure they look sensible. For example it wouldn't make sense to flip the digits horizontally or vertically unless you are expecting to see what in your test data."},{"metadata":{"trusted":true},"cell_type":"code","source":"N_IMAGES = 25\n\nimages = [image for image, label in [train_data[i] for i in range(N_IMAGES)]] \n\nplot_images(images)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The MNIST dataset comes with a training and test set, but not a validation set. We want to use a validation set to check how well our model performs on unseen data. Why don't we just use the test data? We should only be measuring our performance over the test set once, after all training is done. We can think of the validation set as a proxy test set we are allowed to look at as much as we want. \n\nWe create a validation set taking 10% of the training set. **Note:** the validation set should always be created from the training set. Never take the validation set from the test set. When researchers publish research papers they should be comparing performance across the test set and the only way to ensure this is a fair comparison is for all researchers to use the same test set. If the validation set is taken from the test set, then the test set is not the same as everyone else's and the results cannot be compared against each other.\n\nFirst, we have to define the exact number of examples that we want to be in each split of the training/validation sets."},{"metadata":{"trusted":true},"cell_type":"code","source":"VALID_RATIO = 0.9\n\nn_train_examples = int(len(train_data) * VALID_RATIO)\nn_valid_examples = len(train_data) - n_train_examples","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Then, we use the random_split function to take a random 10% of the training set to use as a validation set. The remaining 90% will stay as the training set."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_data, valid_data = data.random_split(train_data, \n                                           [n_train_examples, n_valid_examples])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can print out the number of examples again to check our splits are correct."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Number of training examples: {len(train_data)}')\nprint(f'Number of validation examples: {len(valid_data)}')\nprint(f'Number of testing examples: {len(test_data)}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"One thing to consider is that as the validation set has been created from the training set it has the same transforms as the training set, with the random rotating and cropping. As we want our validation set to act as a proxy for the test set, it should also be fixed, without any random augmentation.\n\nFirst, let's see what 25 of the images within the validation set look like with the training transforms:"},{"metadata":{"trusted":true},"cell_type":"code","source":"N_IMAGES = 25\n\nimages = [image for image, label in [valid_data[i] for i in range(N_IMAGES)]] \n\nplot_images(images)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can now simply replace the validation set's transform by overwriting it with our test transforms from above.\n\nAs the validation set is a Subset of the training set, if we change the transforms of one then by default Torchvision will change the transforms of the other. To stop this from happening we make a deepcopy of the validation data."},{"metadata":{"trusted":true},"cell_type":"code","source":"valid_data = copy.deepcopy(valid_data)\nvalid_data.dataset.transform = test_transforms","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To double check we've correctly replaced the training transforms, we can view the same set of images and notice how they're more central (no random cropping) and have a more standard orientation (no random rotations)."},{"metadata":{"trusted":true},"cell_type":"code","source":"N_IMAGES = 25\n\nimages = [image for image, label in [valid_data[i] for i in range(N_IMAGES)]] \n\nplot_images(images)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we'll define a `DataLoader` for each of the training/validation/test sets. We can iterate over these and they will yield batches of images and labels which we can use to train our model.\n\nWe only need to shuffle our training set as it will be used for stochastic gradient descent and we want the each batch to be different between epochs. As we aren't using the validation or test sets to update our model parameters they do not need to be shuffled.\n\nIdeally, we want to use the biggest batch size that we can. The 64 here is relatively small and can be increased if our hardware can handle it."},{"metadata":{"trusted":true},"cell_type":"code","source":"BATCH_SIZE = 64\n\ntrain_iterator = data.DataLoader(train_data, \n                                 shuffle = True, \n                                 batch_size = BATCH_SIZE)\n\nvalid_iterator = data.DataLoader(valid_data, \n                                 batch_size = BATCH_SIZE)\n\ntest_iterator = data.DataLoader(test_data, \n                                batch_size = BATCH_SIZE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Defining the Model\n\nOur model will be a neural network, specifically a multilayer perceptron (MLP) with two hidden layers. The image below shows the archicture of the model. \n\n![](https://github.com/bentrevett/pytorch-image-classification/blob/master/assets/mlp-mnist.png?raw=1)\n\nSpecifically, first we will flatten our 1x28x28 (1 color channel, 28 pixels height and width) image into a 784 element vector, also called 784 *features*. We flatten our input as MLPs cannot handle two or three dimensional data. Next, the 784 dimensional input is passed through the first hidden layer to transform it into 250 dimensions. Then, another hidden layer, which will transform it to 100 dimensions. Finally, an output layer which will transform it into a 10 dimensional vector. The output dimension should equal the number of classes within your data. Here we have ten digits, 0 - 9, so need our output to be 10 dimensions.\n\nThe transformation between 784 to 250, 250 to 100 and 100 to 10 dimensions are done by `Linear` layers. These are also known as fully connected or affine layers. In these layers every element in one layer is connected to every element in the next. We can think of these elements as *neurons* as this architecture is inspired by how the human brain is made of millions of interconnected nodes, also called neurons. \n\nEach connection between a neuron in one layer and a neuron in the next has a *weight* associated with it. The input to one neuron is sum of the weighted values of all neurons in the previous layer connected to it, plus a weighted bias term, where the bias value is always 1. The neuron then applies an *activation function* to this weighted sum. This activation function is a non-linear function that allows the neural network to learn non-linear functions between inputs and outputs. \n\nWe define our MLP below, which consists of three linear layers. We first take the input batch of images and flatten them so they can be passed into the linear layers. We then pass them through the first linear layer, `input_fc`, which calculates the weighted sum of the inputs, and then apply the *ReLU* (rectified linear unit) activation function elementwise. This result is then passed through another linear layer, `hidden_fc`, again applying the same activation function elementwise. Finally, we pass this through the final linear layer, `output_fc`. We return not only the output but also the second hidden layer as we will do some analysis on it later.\n\nThe ReLU activation function is a popular non-linear function that is simply $max(0, x)$, where $x$ is the weighted sum of the inputs to that neuron. Other activation functions used are hyperbolic tan (tanh) and sigmoid function, however ReLU is the most commonly used.\n\n![](https://github.com/bentrevett/pytorch-image-classification/blob/master/assets/relu.png?raw=1)\n\n\nOne thing to note is that we do not use an activation function on the input directly or on the output. You should never use activation functions directly on the input, i.e. `F.relu(x)`. PyTorch combines activation functions to be applied on the output with the functions which calculate the *loss*, also known as *error* or *cost*, of a neural network. This is done for numerical stability.\n\nWhy did we choose hidden dimensions of 250 and 100 elements? Why did we only have two hidden layers? There is no magic formula to tell us how many layers to use and how many neurons to have in each layer and there is most probably a better set of values. However, the general idea is that neural networks extract features from data. Layers closer to the input learn to extract general features (e.g. lines, curves, edges), whilst later layers combine the features extracted from the previous layer into more high level features (e.g. the intersection of two lines making a cross, multiple curves make a circle). We force our neural network to learn these features by reducing the number of neurons in each layer. This way it has to learn to compress information by extracting only the useful and general features. Thus, we want a neural network with multiple layers and some sort of information compression (reduced number of neurons in subsequent layers)."},{"metadata":{"trusted":true},"cell_type":"code","source":"class MLP(nn.Module):\n    def __init__(self, input_dim, output_dim):\n        super().__init__()\n                \n        self.input_fc = nn.Linear(input_dim, 250)\n        self.hidden_fc = nn.Linear(250, 100)\n        self.output_fc = nn.Linear(100, output_dim)\n        \n    def forward(self, x):\n        \n        #x = [batch size, height, width]\n        \n        batch_size = x.shape[0]\n\n        x = x.view(batch_size, -1)\n        \n        #x = [batch size, height * width]\n        \n        h_1 = F.relu(self.input_fc(x))\n        \n        #h_1 = [batch size, 250]\n\n        h_2 = F.relu(self.hidden_fc(h_1))\n\n        #h_2 = [batch size, 100]\n\n        y_pred = self.output_fc(h_2)\n        \n        #y_pred = [batch size, output dim]\n        \n        return y_pred, h_2","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We'll define our model by creating an instance of it and setting the correct input and output dimensions.\nINPUT_DIM = 28 * 28\nOUTPUT_DIM = 10\n\nmodel = MLP(INPUT_DIM, OUTPUT_DIM)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can also create a small function to calculate the number of trainable parameters (weights and biases) in our model - in case all of our parameters are trainable."},{"metadata":{"trusted":true},"cell_type":"code","source":"def count_parameters(model):\n    return sum(p.numel() for p in model.parameters() if p.requires_grad)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first layer has 784 neurons connected to 250 neurons, so 784*250 weighted connections plus 250 bias terms.\n\nThe second layer has 250 neurons connected to 100 neurons, 250*100 weighted connections plus 100 bias terms.\n\nThe third layer has 100 neurons connected to 10 neurons, 100*10 weighted connections plus 10 bias terms.\n\n$$784 \\cdot 250 + 250 + 250 \\cdot 100 + 100 + 100 \\cdot 10 + 10 = 222,360 $$"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'The model has {count_parameters(model):,} trainable parameters')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training the Model\n\nNext, we'll define our optimizer. This is the algorithm we will use to update the parameters of our model with respect to the loss calculated on the data.\n\nWe aren't going to go into too much detail on how neural networks are trained (see [this](http://neuralnetworksanddeeplearning.com/) article if you want to know how) but the gist is:\n- pass a batch of data through your model\n- calculate the loss of your batch by comparing your model's predictions against the actual labels\n- calculate the gradient of each of your parameters with respect to the loss\n- update each of your parameters by subtracting their gradient multiplied by a small *learning rate* parameter\n\nWe use the *Adam* algorithm with the default parameters to update our model. Improved results could be obtained by searching over different optimizers and learning rates, however default Adam is usually a good starting off point. Check out [this](https://ruder.io/optimizing-gradient-descent/) article if you want to learn more about the different optimization algorithms commonly used for neural networks."},{"metadata":{"trusted":true},"cell_type":"code","source":"optimizer = optim.Adam(model.parameters())","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Loss Function Definition\nThen, we define a *criterion*, PyTorch's name for a loss/cost/error function. This function will take in your model's predictions with the actual labels and then compute the loss/cost/error of your model with its current parameters.\n\n`CrossEntropyLoss` both computes the *softmax* activation function on the supplied predictions as well as the actual loss via *negative log likelihood*. \n\nBriefly, the softmax function is:\n\n$$\\text{softmax }(\\mathbf{x}) = \\frac{e^{x_i}}{\\sum_j e^{x_j}}$$ \n\nThis turns out 10 dimensional output, where each element is an unbounded real number, into a probability distribution over 10 elements. That is, all values are between 0 and 1, and together they all sum to 1. \n\nWhy do we turn things into a probability distribution? So we can use negative log likelihood for our loss function as it expects probabilities. PyTorch calculates negative log likelihood for a single example via:\n\n$$\\text{negative log likelihood }(\\mathbf{\\hat{y}}, y) = -\\log \\big( \\text{softmax}(\\mathbf{\\hat{y}})[y] \\big)$$\n\n$\\mathbf{\\hat{y}}$ is the $\\mathbb{R}^{10}$ output, from our neural network, whereas $y$ is the label, an integer representing the class. The loss is the negative log of the class index of the softmax. For example:\n\n$$\\mathbf{\\hat{y}} = [5,1,1,1,1,1,1,1,1,1]$$\n\n$$\\text{softmax }(\\mathbf{\\hat{y}}) = [0.8585, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157, 0.0157]$$\n\nIf the label was class zero, the loss would be:\n\n$$\\text{negative log likelihood }(\\mathbf{\\hat{y}}, 0) = - \\log(0.8585) = 0.153 \\dots$$\n\nIf the label was class five, the loss would be:\n\n$$\\text{negative log likelihood }(\\mathbf{\\hat{y}}, 5) = - \\log(0.0157) = 4.154 \\dots$$\n\nSo, intuitively, as your model's output corresponding to the correct class index increases your loss decreases."},{"metadata":{"trusted":true},"cell_type":"code","source":"criterion = nn.CrossEntropyLoss()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We then define device. This is used to place your model and data on to a GPU, if you have one."},{"metadata":{"trusted":true},"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We place our model and criterion on to the device by using the .to method."},{"metadata":{"trusted":true},"cell_type":"code","source":"model = model.to(device)\ncriterion = criterion.to(device)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we'll define a function to calculate the accuracy of our model. This takes the index of the highest value for your prediction and compares it against the actual class label. We then divide how many our model got correct by the amount in the batch to calculate accuracy across the batch."},{"metadata":{"trusted":true},"cell_type":"code","source":"def calculate_accuracy(y_pred, y):\n    top_pred = y_pred.argmax(1, keepdim = True)\n    correct = top_pred.eq(y.view_as(top_pred)).sum()\n    acc = correct.float() / y.shape[0]\n    return acc","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We finally define our training loop.\n\nThis will:\n- put our model into `train` mode\n- iterate over our dataloader, returning batches of (image, label)\n- place the batch on to our GPU, if we have one\n- clear the gradients calculated from the last batch\n- pass our batch of images, `x`, through to model to get predictions, `y_pred`\n- calculate the loss between our predictions and the actual labels\n- calculate the accuracy between our predictions and the actual labels\n- calculate the gradients of each parameter\n- update the parameters by taking an optimizer step\n- update our metrics\n\nSome layers act differently when training and evaluating the model that contains them, hence why we must tell our model we are in \"training\" mode. The model we are using here does not use any of those layers, however it is good practice to get used to putting your model in training mode."},{"metadata":{"trusted":true},"cell_type":"code","source":"def train(model, iterator, optimizer, criterion, device):\n    \n    epoch_loss = 0\n    epoch_acc = 0\n    \n    model.train()\n    \n    for (x, y) in iterator:\n        \n        x = x.to(device)\n        y = y.to(device)\n        \n        optimizer.zero_grad()\n                \n        y_pred, _ = model(x)\n        \n        loss = criterion(y_pred, y)\n        \n        acc = calculate_accuracy(y_pred, y)\n        \n        loss.backward()\n        \n        optimizer.step()\n        \n        epoch_loss += loss.item()\n        epoch_acc += acc.item()\n        \n    return epoch_loss / len(iterator), epoch_acc / len(iterator)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The evaluation loop is similar to the training loop. The differences are:\n- we put our model into evaluation mode with `model.eval()`\n- we wrap the iterations inside a `with torch.no_grad()`\n- we do not zero gradients as we are not calculating any\n- we do not calculate gradients as we are not updating parameters\n- we do not take an optimizer step as we are not calculating gradients\n\n`torch.no_grad()` ensures that gradients are not calculated for whatever is inside the `with` block. As our model will not have to calculate gradients it will be faster and use less memory. "},{"metadata":{"trusted":true},"cell_type":"code","source":"def evaluate(model, iterator, criterion, device):\n    \n    epoch_loss = 0\n    epoch_acc = 0\n    \n    model.eval()\n    \n    with torch.no_grad():\n        \n        for (x, y) in iterator:\n\n            x = x.to(device)\n            y = y.to(device)\n\n            y_pred, _ = model(x)\n\n            loss = criterion(y_pred, y)\n\n            acc = calculate_accuracy(y_pred, y)\n\n            epoch_loss += loss.item()\n            epoch_acc += acc.item()\n        \n    return epoch_loss / len(iterator), epoch_acc / len(iterator)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# The final step before training is to define a small function to tell us how long an epoch took.\ndef epoch_time(start_time, end_time):\n    elapsed_time = end_time - start_time\n    elapsed_mins = int(elapsed_time / 60)\n    elapsed_secs = int(elapsed_time - (elapsed_mins * 60))\n    return elapsed_mins, elapsed_secs","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We're finally ready to train!\n\nDuring each epoch we calculate the training loss and accuracy, followed by the validation loss and accuracy. We then check if the validation loss achieved is the best validation loss we have seen. If so, we save our model's parameters (called a state_dict)."},{"metadata":{"trusted":true},"cell_type":"code","source":"EPOCHS = 10\n\nbest_valid_loss = float('inf')\n\nfor epoch in range(EPOCHS):\n    \n    start_time = time.time()\n    \n    train_loss, train_acc = train(model, train_iterator, optimizer, criterion, device)\n    valid_loss, valid_acc = evaluate(model, valid_iterator, criterion, device)\n    \n    if valid_loss < best_valid_loss:\n        best_valid_loss = valid_loss\n        torch.save(model.state_dict(), 'tut1-model.pt')\n    \n    end_time = time.time()\n\n    epoch_mins, epoch_secs = epoch_time(start_time, end_time)\n    \n    print(f'Epoch: {epoch+1:02} | Epoch Time: {epoch_mins}m {epoch_secs}s')\n    print(f'\\tTrain Loss: {train_loss:.3f} | Train Acc: {train_acc*100:.2f}%')\n    print(f'\\t Val. Loss: {valid_loss:.3f} |  Val. Acc: {valid_acc*100:.2f}%')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Afterwards, we load our the parameters of the model that achieved the best validation loss and then use this to evaluate our model on the test set."},{"metadata":{"trusted":true},"cell_type":"code","source":"model.load_state_dict(torch.load('tut1-model.pt'))\n\ntest_loss, test_acc = evaluate(model, test_iterator, criterion, device)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Test Loss: {test_loss:.3f} | Test Acc: {test_acc*100:.2f}%')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our model achieves 98% accuracy on the test set.\n\nThis can be improved by tweaking hyperparameters, e.g. number of layers, number of neurons per layer, optimization algorithm used, learning rate, etc."},{"metadata":{},"cell_type":"markdown","source":"# Examining the Model\n\nNow we've trained our model there's a few things we can look at. Most of these are simple exploratory analysis, but they can offer some insights into your model.\n\nAn important thing to do is check what examples your model gets wrong and ensure that they're reasonable mistakes.\n\nThe function below will return the model's predictions over a given dataset. It will return the inputs (image) the outputs (model predictions) and the ground truth labels."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_predictions(model, iterator, device):\n\n    model.eval()\n\n    images = []\n    labels = []\n    probs = []\n\n    with torch.no_grad():\n\n        for (x, y) in iterator:\n\n            x = x.to(device)\n\n            y_pred, _ = model(x)\n\n            y_prob = F.softmax(y_pred, dim = -1)\n            top_pred = y_prob.argmax(1, keepdim = True)\n\n            images.append(x.cpu())\n            labels.append(y.cpu())\n            probs.append(y_prob.cpu())\n\n    images = torch.cat(images, dim = 0)\n    labels = torch.cat(labels, dim = 0)\n    probs = torch.cat(probs, dim = 0)\n\n    return images, labels, probs","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We can then get these predictions and, by taking the index of the highest predicted probability, get the predicted labels.\nimages, labels, probs = get_predictions(model, test_iterator, device)\n\npred_labels = torch.argmax(probs, 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Next, for each of our examples, we can check if our predicted label matches our actual label.\ncorrects = torch.eq(labels, pred_labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can then loop through all of the examples over our model's predictions and store all the examples the model got incorrect into an array.\n\nThen, we sort these incorrect examples by how confident they were, with the most confident being first."},{"metadata":{"trusted":true},"cell_type":"code","source":"incorrect_examples = []\n\nfor image, label, prob, correct in zip(images, labels, probs, corrects):\n    if not correct:\n        incorrect_examples.append((image, label, prob))\n\nincorrect_examples.sort(reverse = True, key = lambda x: torch.max(x[2], dim = 0).values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can then plot the incorrectly predicted images along with how confident they were on the actual label and how confident they were at the incorrect label."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_most_incorrect(incorrect, n_images):\n\n    rows = int(np.sqrt(n_images))\n    cols = int(np.sqrt(n_images))\n\n    fig = plt.figure(figsize = (20, 10))\n    for i in range(rows*cols):\n        ax = fig.add_subplot(rows, cols, i+1)\n        image, true_label, probs = incorrect[i]\n        true_prob = probs[true_label]\n        incorrect_prob, incorrect_label = torch.max(probs, dim = 0)\n        ax.imshow(image.view(28, 28).cpu().numpy(), cmap='bone')\n        ax.set_title(f'true label: {true_label} ({true_prob:.3f})\\n' \\\n                     f'pred label: {incorrect_label} ({incorrect_prob:.3f})')\n        ax.axis('off')\n    fig.subplots_adjust(hspace= 0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Below we can see the 25 images the model got incorrect and was most confident about.\n\nA lot of these digits are irregular so it's difficult for the model to do well on these. The images that do look fine, if you squint you can sort of see why the model got it wrong.\n\nWhy is the neural network so confident on the irregular digits? Surely if it's a weird looking digit then the output of the softmax should be close to evenly distributed across a few digits the model isn't sure about, right? Well, no. The model has been trained to only be incredibly confident about its predictions and thus when it sees an image it will always be confident about what it is."},{"metadata":{"trusted":true},"cell_type":"code","source":"N_IMAGES = 25\n\nplot_most_incorrect(incorrect_examples, N_IMAGES)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Another thing we can do is get the output and intermediate representations from the model and try to visualize them.\n\nThe function below loops through the provided dataset and gets the output from the model and the intermediate representation from the layer before that, the second hidden layer."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_representations(model, iterator, device):\n\n    model.eval()\n\n    outputs = []\n    intermediates = []\n    labels = []\n\n    with torch.no_grad():\n        \n        for (x, y) in iterator:\n\n            x = x.to(device)\n\n            y_pred, h = model(x)\n\n            outputs.append(y_pred.cpu())\n            intermediates.append(h.cpu())\n            labels.append(y)\n        \n    outputs = torch.cat(outputs, dim = 0)\n    intermediates = torch.cat(intermediates, dim = 0)\n    labels = torch.cat(labels, dim = 0)\n\n    return outputs, intermediates, labels","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We run the function to get the representations."},{"metadata":{"trusted":true},"cell_type":"code","source":"outputs, intermediates, labels = get_representations(model, train_iterator, device)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PCA"},{"metadata":{},"cell_type":"markdown","source":"The data we want to visualize is in ten dimensions and 100 dimensions. We want to get this down to two dimensions so we can actually plot it.\n\nThe first technique we'll use is PCA (principal component analysis). First, we'll define a function to calculate the PCA of our data and then we'll define a function to plot it."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_pca(data, n_components = 2):\n    pca = decomposition.PCA()\n    pca.n_components = n_components\n    pca_data = pca.fit_transform(data)\n    return pca_data","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_representations(data, labels, n_images = None):\n    if n_images is not None:\n        data = data[:n_images]\n        labels = labels[:n_images]\n    fig = plt.figure(figsize = (10, 10))\n    ax = fig.add_subplot(111)\n    scatter = ax.scatter(data[:, 0], data[:, 1], c = labels, cmap = 'tab10')\n    handles, labels = scatter.legend_elements()\n    legend = ax.legend(handles = handles, labels = labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First, we plot the representations from the ten dimensional output layer, reduced down to two dimensions."},{"metadata":{"trusted":true},"cell_type":"code","source":"output_pca_data = get_pca(outputs)\nplot_representations(output_pca_data, labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Next, we'll plot the outputs of the second hidden layer.\n\nThe clusters seem similar to the one above. In fact if we rotated the below image anti-clockwise it wouldn't be too far off the PCA of the output representations."},{"metadata":{"trusted":true},"cell_type":"code","source":"intermediate_pca_data = get_pca(intermediates)\nplot_representations(intermediate_pca_data, labels)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# t-SNE\nAn alternative to PCA is t-SNE (t-distributed stochastic neighbor embedding). \n\nThis is commonly thought of as being \"better\" than PCA, although it can be [misinterpretend](https://distill.pub/2016/misread-tsne/)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_tsne(data, n_components = 2, n_images = None):\n    if n_images is not None:\n        data = data[:n_images]\n    tsne = manifold.TSNE(n_components = n_components, random_state = 0)\n    tsne_data = tsne.fit_transform(data)\n    return tsne_data","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"t-SNE is very slow, so we only compute it on a subset of the representations.\n\nThe classes look very well separated, and it is possible to use k-NN on this representation to achieve decent accuracy."},{"metadata":{"trusted":true},"cell_type":"code","source":"N_IMAGES = 5_000\n\noutput_tsne_data = get_tsne(outputs, n_images = N_IMAGES)\nplot_representations(output_tsne_data, labels, n_images = N_IMAGES)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We plot the intermediate representations on the same subset.\n\nAgain, the classes look well separated, though less so than the output representations. This is because these representations are intermediate features that the neural network has extracted and will use them in further layers to weigh up the evidence of what digit is in the image. Hence, in theory, the classes should become more separated the closer we are to the output layer, which is exactly what we see here."},{"metadata":{"trusted":true},"cell_type":"code","source":"intermediate_tsne_data = get_tsne(intermediates, n_images = N_IMAGES)\nplot_representations(intermediate_tsne_data, labels, n_images = N_IMAGES)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Another experiment we can do is try and generate fake digits.\n\nThe function below will repeatedly generate random noise and feed it through the model and find the most confidently generated digit for the desired class."},{"metadata":{"trusted":true},"cell_type":"code","source":"def imagine_digit(model, digit, device, n_iterations = 50_000):\n\n    model.eval()\n    \n    best_prob = 0\n    best_image = None\n\n    with torch.no_grad():\n    \n        for _ in range(n_iterations):\n\n            x = torch.randn(32, 28, 28).to(device)\n\n            y_pred, _ = model(x)\n\n            preds = F.softmax(y_pred, dim = -1)\n\n            _best_prob, index = torch.max(preds[:,digit], dim = 0)\n\n            if _best_prob > best_prob:\n                best_prob = _best_prob\n                best_image = x[index]\n\n    return best_image, best_prob","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DIGIT = 3\n\nbest_image, best_prob = imagine_digit(model, DIGIT, device)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looking at the best probability achieved, we have a digit that the model is 100% confident is a three."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(f'Best image probability: {best_prob.item()*100:.2f}%')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Unfortunately the imagined perfect three just looks like random noise.\n\nAs mentioned before, the model has only been trained to be incredibly confident with its predictions so when faced with random noise will try and classify it as something.\n\nIt is also possible that the model is overfitting on the training set - that there is a common pattern in handwritten threes within the training set but it's not the pattern we want our model to learn."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(best_image.cpu().numpy(), cmap = 'bone')\nplt.axis('off');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally, we can plot the weights in the first layer of our model.\n\nThe hope is that there's maybe one neuron in this first layer that's learned to look for certain patterns in the input and thus has high weight values indicating this pattern. If we then plot these weights we should see these patterns."},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_weights(weights, n_weights):\n\n    rows = int(np.sqrt(n_weights))\n    cols = int(np.sqrt(n_weights))\n\n    fig = plt.figure(figsize = (20, 10))\n    for i in range(rows*cols):\n        ax = fig.add_subplot(rows, cols, i+1)\n        ax.imshow(weights[i].view(28, 28).cpu().numpy(), cmap = 'bone')\n        ax.axis('off')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Looking at these weights we see a few of them look like random noise but some of them do have weird patterns within them. These patterns show \"ghostly\" image looking shapes, but are clearly not images."},{"metadata":{"trusted":true},"cell_type":"code","source":"N_WEIGHTS = 25\n\nweights = model.input_fc.weight.data\n\nplot_weights(weights, N_WEIGHTS)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Conclusions\n\nIn this notebook we have shown: \n- loading Torchvision datasets\n- loading transforms to augment and normalize our data\n- defining a MLP\n- training a model to achieve >97% accuracy\n- viewing our model's mistakes\n- visualizing our data in lower dimensions with PCA and t-SNE\n- generating fake digits\n- viewing the learned weights of our model"},{"metadata":{},"cell_type":"markdown","source":"# <font color='Red'> Thankyou For Reading And Do Upvote If You Liked the Notebook!!!!"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}