final_project_2022-2023.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final project, Numerical Analysis 2022-2023\n",
    "\n",
    "\n",
    "## Project description\n",
    "\n",
    "In this project, we would like to compare the performance of some embarassingly simple algorithms to solve a classification problem based on the MNIST database. \n",
    "\n",
    "The abstract aim of the program is to write a function:\n",
    "\n",
    "```\n",
    "result = classify(image)\n",
    "```\n",
    "\n",
    "that takes as input a small grey scale image of a hand-written digit (from the MNIST database), and returns the digit corresponding to the content of the image.\n",
    "\n",
    "An example of the images we'll be working on is the following:\n",
    "\n",
    "![mnist examples](https://m-alcu.github.io/assets/mnist.png)\n",
    "\n",
    "Some background on the MNIST database (from wikipedia):\n",
    "\n",
    "\n",
    "## MNIST database\n",
    "\n",
    "*From Wikipedia, the free encyclopedia*\n",
    "\n",
    "The MNIST database (Modified National Institute of Standards and Technology database) is a large database of handwritten digits that is commonly used for training various image processing systems. The database is also widely used for training and testing in the field of machine learning. It was created by \"re-mixing\" the samples from NIST's original datasets. The creators felt that since NIST's training dataset was taken from American Census Bureau employees, while the testing dataset was taken from American high school students, it was not well-suited for machine learning experiments. Furthermore, the black and white images from NIST were normalized to fit into a 28x28 pixel bounding box and anti-aliased, which introduced grayscale levels.\n",
    "\n",
    "## MNIST sample images.\n",
    "\n",
    "The MNIST database contains 60,000 training images and 10,000 testing images. Half of the training set and half of the test set were taken from NIST's training dataset, while the other half of the training set and the other half of the test set were taken from NIST's testing dataset. There have been a number of scientific papers on attempts to achieve the lowest error rate; one paper, using a hierarchical system of convolutional neural networks, manages to get an error rate on the MNIST database of 0.23%. The original creators of the database keep a list of some of the methods tested on it. In their original paper, they use a support vector machine to get an error rate of 0.8%. An extended dataset similar to MNIST called EMNIST has been published in 2017, which contains 240,000 training images, and 40,000 testing images of handwritten digits and characters.\n",
    "\n",
    "## Algorithm\n",
    "\n",
    "We start by defining the distance between two images. Ideally, a distance function between two images is zero when the images are the same, and greater than zero when the images are different. \n",
    "\n",
    "The bigger the distance, the more different the images should be. Ideally, the distance between an image of the number `9` should be closer to an image of the number `8` than to an image of the number `1` (the digits `9` and `8`, as images, differ by the fact that the first has one closed loop, while the second has two closed loops, while the digit `1` is mostly a straight line). Two different images representing the same number should be even closer (i.e., the distance function should return a \"small\" number).\n",
    "\n",
    "Given a distance and a training set of images for which we know everything, the simplest algorithm we can think of to classify an image `z`, is the following: given a set of train images (`x_train`) for which we know the digit they represent (`y_train`), measure the distance between `z` and all images in `x_train`, and classify the image `z` to represent the same digit of the image that is closest to `z` in `x_train`:\n",
    "\n",
    "Parameters of the algorithm:\n",
    "\n",
    "- `x_train`\n",
    "- `y_train`\n",
    "- a distance function `dist`  \n",
    "\n",
    "Input of the function\n",
    "\n",
    "- `z`\n",
    "\n",
    "Output of the function\n",
    "\n",
    "- `digit`\n",
    "\n",
    "where \n",
    "\n",
    "```\n",
    "def classify(z):\n",
    "    all_distances = array([dist(x, z) for x in x_train])\n",
    "    digit = y_train[argmin(all_distances)]\n",
    "    return digit\n",
    "```\n",
    "\n",
    "We will experiment with different distances, and we will try to improve the algorithm above in a step by step fashon.\n",
    "\n",
    "## Data description\n",
    "\n",
    "Each image in the MNIST dataset represents a hand written digit, in the form of a matrix of `28x28` values between zero and one, representing gray scale values (zero = white, one = black).\n",
    "\n",
    "We use an array of `60.000x28x28` floating point values to collect all training images, and an array of `60.000` digits containing the (correct) value of the training digits (between 0 and 9 inclusive).\n",
    "\n",
    "The testing images are instead collected into two arrays of size `10.000x28x28` and `10.0000` respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "\n",
    "# # input image dimensions\n",
    "# img_rows, img_cols = 28, 28\n",
    "\n",
    "# # Uncomment the following lines if you have keras installed. Otherwise you can \n",
    "# # use the file I uploaded: mnist.npz\n",
    "# import keras\n",
    "# from keras.datasets import mnist\n",
    "\n",
    "# # the data, split between train and test sets\n",
    "# (x_train, y_train), (x_test, y_test) = mnist.load_data()\n",
    "\n",
    "# if K.image_data_format() == 'channels_first':\n",
    "#     x_train = x_train.reshape(x_train.shape[0], 1, img_rows, img_cols)\n",
    "#     x_test = x_test.reshape(x_test.shape[0], 1, img_rows, img_cols)\n",
    "#     input_shape = (img_rows, img_cols)\n",
    "# else:\n",
    "#     x_train = x_train.reshape(x_train.shape[0], img_rows, img_cols)\n",
    "#     x_test = x_test.reshape(x_test.shape[0], img_rows, img_cols)\n",
    "#     input_shape = (img_rows, img_cols)\n",
    "\n",
    "# x_train = x_train.astype('float32')\n",
    "# x_test = x_test.astype('float32')\n",
    "# x_train /= 255\n",
    "# x_test /= 255\n",
    "\n",
    "# np.savez_compressed('mnist.npz', x_train, y_train, x_test, y_test)\n",
    "\n",
    "arc = load('mnist.npz')\n",
    "\n",
    "x_train = arc['arr_0']\n",
    "y_train = arc['arr_1']\n",
    "x_test  = arc['arr_2']\n",
    "y_test  = arc['arr_3']\n",
    "\n",
    "print(x_train.shape, y_train.shape)\n",
    "print(x_test.shape, y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting one image\n",
    "\n",
    "How do we plot the images? `pyplot`, which has been imported by the first line of the previous cell, contains a command called `imshow`, that can be used to plot images. \n",
    "\n",
    "In this case we know it is a greyscale image, with zero representing white and one representing black, so we use a colormap that goes from white to black, i.e., `gray_r` where `_r` stands for reversed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show image number 15, and write in the title what digit it should correspond to\n",
    "N=15\n",
    "imshow(x_train[N], cmap='gray_r')\n",
    "_ = title('Hand written digit '+str(y_train[N]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**IF YOU DON'T HAVE ENOUGH COMPUTATIONAL POWER, RUN THE EXERCISES ONLY UP TO WHAT IS SUSTAINABLE FOR YOUR PC**\n",
    "\n",
    "General guidelines:\n",
    "\n",
    "- Time all functions you construct, and try to make them run as fast as possible by precomputing anything that can be precomputed\n",
    "- Extra points are gained if you reduce the complexity of the given algorithms in any possible way, for example by exploiting linearity, etc.\n",
    "- If something takes too long to execute, make sure you time it on a smaller set of input data, and give estimates of how long it would take to run the full thing (without actually running it). Plot only the results you manage to run on your PC.\n",
    "\n",
    "# Assignment 1\n",
    "\n",
    "Implement the following distance functions\n",
    "\n",
    "- d_infty $$ d_{\\infty}(a,b) := \\max_{i,j} |b_{ij}-a_{ij}|$$\n",
    "- d_one $$ d_1(a,b) := \\sum_{i,j} |b_{ij}-a_{ij}|$$\n",
    "- d_two $$ d_2(a,b) := \\sqrt{\\sum_{i,j} |b_{ij}-a_{ij}|^2}$$\n",
    "\n",
    "that take two `(28,28)` images in input, and return a non-negative number.\n",
    "\n",
    "# Assignment 2\n",
    "\n",
    "Write a function that, given a number `N`, and a distance function `dist`, computes the distance matrix D of shape `(N,N)` between the first `N` entries of `x_train`:\n",
    "\n",
    "```\n",
    "D[i,j] = dist(x_train[i], x_train[j])\n",
    "```\n",
    "\n",
    "performing the **minimum** number of operations (i.e., avoid computing a distance if it has already been computed before, i.e., keep in mind that dist(a,b) = dist(b,a)).\n",
    "\n",
    "# Assignment 3\n",
    "\n",
    "Compute and plot the three distance matrices\n",
    "\n",
    "- Dinfty\n",
    "- D1\n",
    "- D2\n",
    "\n",
    "for the first 100 images of the training set, using the function `imshow` applied to the three matrices\n",
    "\n",
    "# Assignment 4\n",
    "\n",
    "Using only a distance matrix, apply the algorithm described above and compute the efficency of the algorithm, i.e., write a function that:\n",
    "\n",
    "Given a distance matrix with shape `(N,N)`, constructed on the first `N` samples of the `x_train` set, count the number of failures of the **leave one out** strategy, i.e., \n",
    "\n",
    "- set `error_counter` to zero\n",
    "\n",
    "- for every line `i` of the matrix:\n",
    "\n",
    "    - find the index `j` (different from `i`) for which `D[i,k] >= D[i,j]` for all `k` different from `i` and `j`.\n",
    "\n",
    "    - if `y_train[j]` is different from `y_train[i]`, increment by one `error_counter`\n",
    "\n",
    "- return the error: error_counter/N\n",
    "\n",
    "- apply the function above to the 3 different distance matrices you computed before\n",
    "\n",
    "\n",
    "# Assignment 5\n",
    "\n",
    "Run the algorithm implemented above for N=100,200,400,800,1600 on the three different distances, and plot the three error rate as a function of N (i.e., compute the distance matrix, and compute the efficiency associated to the distance matrix).\n",
    "\n",
    "You should get an error like:\n",
    "```\n",
    "[[ 0.58      0.17      0.17    ]\n",
    " [ 0.52      0.145     0.135   ]\n",
    " [ 0.4425    0.15      0.135   ]\n",
    " [ 0.4       0.145     0.12875 ]\n",
    " [ 0.369375  0.1025    0.09375 ]]\n",
    "```\n",
    "where each column represents a different norm.\n",
    "\n",
    "** In the next assignments, optional points are given if you manage to make the algorithm run faster, by pre-computing everything you can precompute in advance **\n",
    "\n",
    "# Assignment 6\n",
    "\n",
    "In principle, it should be possible to decrease the error by using a better norm. From the table above, it is clear that the L2 distance works better than the L1 distance, which works better than the Linfty distance.\n",
    "\n",
    "However, *none of these distances exploit the fact that the image is a two-dimensional object*, and that there is information also in the **neighboring** information of the pixels. \n",
    "\n",
    "One way to exploit this, is to interpret the image as a continuous function with values between zero and one, defined on a square domain `\\Omega=[0,27]x[0,27]`.\n",
    "\n",
    "$$ f: \\Omega \\to R $$\n",
    "\n",
    "- Implement a function that computes an approximation of the $H^1$ norm distance on the renormalized images. Given two images $f_1$ and $f_2$\n",
    "    - Compute $$a = \\frac{f_1}{\\int_\\Omega f_1}$$, $$b=\\frac{f_2}{\\int_\\Omega f_2}$$\n",
    "    - Define the $H^1$ distance as\n",
    "    $$\n",
    "    d_{H^1}(f_1,f_2) := \\sqrt{\\int_\\Omega |\\nabla(a-b)|^2+ (a-b)^2}\n",
    "    $$\n",
    "    using the algorithm you prefer (or the library you prefer) to compute the gradients and the integrals. Notice that $\\nabla f = (\\partial f/\\partial x, \\partial f/\\partial y)$ is a vector valued function, and $|\\nabla g|^2 := (\\partial g/\\partial x)^2 + (\\partial g/\\partial y)^2$\n",
    "\n",
    "- Compute the distance matrix and the efficiency for this distance for N=100,200,400,800,1600\n",
    "\n",
    "## Assigment 7\n",
    "\n",
    "Use the `BallTree` algorithm (https://en.wikipedia.org/wiki/Ball_tree), from the `sklearn` package, and construct a tree data structure **that uses a distance from the steps defined above**.\n",
    "\n",
    "For N in `[3200, 6400]` and for each distance defined above\n",
    "\n",
    "- Build a tree using the first N entries of the training set `x_train`\n",
    "- Construct a function that tests the efficiency on all the entries of the test set `x_test`:\n",
    "    - for any image in `x_test`, call it `x_test[i]`, query the tree for the nearest neighbor (call it `k`), and assign as predicted digit the digit of the `x_train[k]` image, i.e., `y_train[k]`\n",
    "    - check if `y_train[k]` is equal to the corresponding entry in `y_test[i]`. If not, increment a counter of the error\n",
    "    - return the efficiency, i.e., `error_counter/len(x_test)`\n",
    "- Plot, in a single graph, the error of each distance as a function of `N` (including labels, titles, etc.)\n",
    "\n",
    "- Once you have the tree, experiment with different nearest neighbor algorithms, i.e., instead of taking only one nearest neighbor, take a larger number (a small number of your choice), and instead of returning the single closest digit, return the one with the largest number of occurrences. Plot the same graph you gave before, and see if you gain an improvement. Motivate all choices you have to make to get to the final answer.\n",
    "\n",
    "\n",
    "**IF YOU DON'T HAVE ENOUGH COMPUTATIONAL POWER, RUN THE EXERCISES ONLY UP TO WHAT IS SUSTAINABLE FOR YOUR PC**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}