one-fell-swoop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dataclasses import dataclass # for neural network construction\n",
    "import pickle                     # to import MNIST data\n",
    "import gzip                       # to import MNIST data\n",
    "import random                     # to initialize weights and biases\n",
    "import numpy as np                # for all needed math\n",
    "from PIL import Image, ImageOps   # for image file processing\n",
    "from time import time             # for performance measurement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic math functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sigmoid function \n",
    "\n",
    "It is used as the activation function in all layers of our network and is defined as:\n",
    "\n",
    "$$\\sigma(z) = \\frac{1}{1 + e^{-z}}$$\n",
    "\n",
    "where $z$ is an `ndarray` (vector) and it will apply the formula element-wise and return another `ndarray` of the same shape as $z$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(z):\n",
    "    return 1.0 / (1.0 + np.exp(-z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivative of the sigmoid function\n",
    "\n",
    "Is used in the backpropagation algorithm and is defined as:\n",
    "\n",
    "$$ \\sigma^\\prime(z) = \\left(\\frac{1}{1 + e^{-z}}\\right)^\\prime = \\left((1 + e^{-z})^{-1}\\right)^\\prime $$\n",
    "\n",
    "using the _power rule_ we can rewrite this as:\n",
    "\n",
    "$$ -\\frac{(1 + e^{-z})^\\prime}{(1 + e^{-z})^{2}} $$\n",
    "\n",
    "then using the _sum rule_ rewrite it as (taking the derivative of the nominator):\n",
    "\n",
    "$$ -(1 + e^{-z})^{-2}(e^{-z})^\\prime$$\n",
    "\n",
    "applying the _exponential rule_ we get:\n",
    "\n",
    "$$ -(1 + e^{-z})^{-2}e^{-z}(-z)^\\prime $$ \n",
    "\n",
    "derivative of $-z$ is simply $-1$ that cancels the minus sign, so we now have:\n",
    "\n",
    "$$ (1 + e^{-z})^{-2}e^{-z}$$\n",
    "\n",
    "which is equivalent to:\n",
    "\n",
    "$$ \\frac{e^{-z}}{(1 + e^{-z})^{2}} $$ \n",
    "\n",
    "that can be rewritten as:\n",
    "\n",
    "$$ \\frac{1}{1 + e^{-z}}\\frac{e^{-z}}{1 + e^{-z}}$$\n",
    "\n",
    "and now with a sleight of hand (by adding and subtracting $1$) we can make it to be:\n",
    "\n",
    "$$ \\frac{1}{1 + e^{-z}}\\frac{e^{-z} + 1 - 1}{1 + e^{-z}} = \\frac{1}{1 + e^{-z}}\\left(\\frac{1 + e^{-z}}{1 + e^{-z}} - \\frac{1}{1 + e^{-z}}\\right)$$\n",
    "\n",
    "and after simplifying it, we end up with the formula:\n",
    "\n",
    "$$\\frac{1}{1 + e^{-z}}\\left(1 - \\frac{1}{1 + e^{-z}}\\right) = \\sigma(z)(1 - \\sigma(z)) $$\n",
    "\n",
    "that the code below implements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid_prime(z):\n",
    "    return sigmoid(z) * (1 - sigmoid(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivative of the cost function\n",
    "\n",
    "We use simple mean square error as our cost function:\n",
    "\n",
    "$$ C(a_L) = \\frac{1}{2}\\|y(x)-a_L(x)\\|^2 $$\n",
    "\n",
    "where the cost function is for a _single_ input $x$, and $a_L$ is the activation function (i.e. output) of the last layer of the network.\n",
    "\n",
    "And its derivative (which is used in the backpropagation algorithm) is defined as:\n",
    "\n",
    "$$ \\left[\\frac{dC_x}{da_L}\\right] \\equiv C^\\prime(a_L)$$\n",
    "\n",
    "We can rewrite the cost function as:\n",
    "\n",
    "$$ C(a_L) = \\frac{1}{2}(y^2 - 2ya_L + a_L^2) $$\n",
    "\n",
    "and therefore its derivative is:\n",
    "\n",
    "$$ C^\\prime(a_L) = \\frac{1}{2}(2a_L - 2y) = a_L - y $$\n",
    "\n",
    "Note that again $x, y,$ and $a_L$ are `ndarray`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost_derivative(output_activations, y):\n",
    "    return (output_activations - y) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The neural network\n",
    "\n",
    "The network is not implemented as a full-blown class, as it was in the original code. It is done this way to allow markdown annotations of each function (because there is no simple or elegant way to split the class implementation into multiple cells of the notebook). It is initialized with a list of numbers, each number representing the number of neurons in the corresponding layer of the network. \n",
    "\n",
    "`init_network` returns a _dataclass_ with three elements: num_layers, biases, and weights. The latter two are lists of `num_layers-1` elements of the `ndarray` type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class Network:\n",
    "    num_layers: int\n",
    "    biases: list\n",
    "    weights: list\n",
    "\n",
    "def init_network(layers):\n",
    "    \n",
    "    return Network(\n",
    "        len(layers),\n",
    "        \n",
    "        # input layer doesn't have biases\n",
    "        [np.random.randn(y, 1) for y in layers[1:]],\n",
    "        \n",
    "        # there are no (weighted) connections into input layer or out of the output layer\n",
    "        [np.random.randn(y, x) for x, y in zip(layers[:-1], layers[1:])]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forward propagation function\n",
    "\n",
    "Calculates activation vector for each layer and returns the last activation vector as `ndarray`.\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `nn`: the neural network object returned by the `init_network` function.\n",
    "- `a`: the network input as `ndarray`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feedforward(nn, a):\n",
    "    for b, w in zip(nn.biases, nn.weights):\n",
    "        a = sigmoid(np.dot(w, a) + b)\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation function \n",
    "\n",
    "Returns the number of test inputs for which the neural network outputs the correct result. Note that the neural    network's output is assumed to be the index of whichever neuron in the final layer has the highest activation (one-hot encoding).\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `nn`: the neural network object returned by the `init_network` function.\n",
    "- `test_data`: a list of tuples (x, y) where x is the input and y is the correct output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate(nn, test_data):\n",
    "    test_results = [(np.argmax(feedforward(nn, x)), y) for (x, y) in test_data]\n",
    "    \n",
    "    return sum(int(x == y) for (x, y) in test_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning function\n",
    "\n",
    "This is the main function, however the actual work is done by `stochastic_gradient_descent`, which this function invokes for each mini-batch of the training data.\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `nn`: the neural network object returned by the `init_network` function.\n",
    "- `training_data`: a list of tuples (x, y) representing the training inputs and (known) correct outputs.\n",
    "- `epochs`: number of epochs of training. During each epoch we go through all training data, but in a different order.\n",
    "- `mini_batch_size`: the number of samples in a mini-batch.\n",
    "- `learning_rate`: a hyper-parameter, 0.3 is the recommended value.\n",
    "- `test_data`: another list of tuples (x, y) representing the test inputs and known outputs.\n",
    "\n",
    "If `test_data` is provided then the network will be evaluated against the test data after each epoch, and partial progress printed out. This is useful for tracking progress, but slows things down."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def learn(nn, training_data, epochs, mini_batch_size, learning_rate, test_data = None):\n",
    "    n = len(training_data)\n",
    "\n",
    "    for j in range(epochs):\n",
    "        random.shuffle(training_data) # that's where \"stochastic\" comes from\n",
    "\n",
    "        mini_batches = [\n",
    "            training_data[k: k + mini_batch_size] for k in range(0, n, mini_batch_size)\n",
    "        ]\n",
    "\n",
    "        for mini_batch in mini_batches:\n",
    "            batch_stochastic_gradient_descent(nn, mini_batch, learning_rate) # that's where learning really happes\n",
    "\n",
    "        if test_data:\n",
    "            print('Epoch {0}: accuracy {1}%'.format(f'{j + 1:2}', 100.0 * evaluate(nn, test_data) / len(test_data)))\n",
    "        else:\n",
    "            print('Epoch {0} complete.'.format(f'{j + 1:2}'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Worker function\n",
    "\n",
    "\n",
    "Updates the network's weights and biases by applying stochastic gradient descent using backpropagation to \n",
    "a single mini-batch according to these equations:\n",
    "\n",
    "$$w = w^\\prime - \\frac{\\eta}{m}\\nabla{w}$$\n",
    "\n",
    "$$b = b^\\prime - \\frac{\\eta}{m}\\nabla{b}$$\n",
    "\n",
    "where $m$ is the number of samples in the mini-batch, and $\\eta$ is the learning rate.\n",
    "\n",
    "$\\nabla{w}$ and $\\nabla{b}$ are calculated by consequently adding $\\delta\\nabla{w}$ and $\\delta\\nabla{b}$ for each sample in the mini-batch, and $\\delta\\nabla{w}$ and $\\delta\\nabla{b}$ are computed by the backpropagation function.\n",
    "\n",
    "Parameters:\n",
    "- `nn`: the neural network object returned by the `init_network` function.\n",
    "- `mini_batch`: a list of training data tuples (x, y) representing the training inputs and (known) correct outputs.\n",
    "- `eta`: the learning rate hyper-parameter $\\eta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_stochastic_gradient_descent(nn, mini_batch, eta):\n",
    "    # \"nabla\" is the gradient symbol\n",
    "    nabla_b = [np.zeros(b.shape) for b in nn.biases]\n",
    "    nabla_w = [np.zeros(w.shape) for w in nn.weights]\n",
    "    \n",
    "    # process the whole mini-batch in one fell swoop instead of iterating through each (x, y) tuple\n",
    "    nabla_b, nabla_w = batch_backprop(nn, mini_batch)\n",
    "        \n",
    "    nn.weights = [w - (eta / len(mini_batch)) * nw for w, nw in zip(nn.weights, nabla_w)]\n",
    "    nn.biases  = [b - (eta / len(mini_batch)) * nb for b, nb in zip(nn.biases, nabla_b)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backpropagation\n",
    "\n",
    "This function operates on the whole mini-batch at once.\n",
    "\n",
    "Backpropagation function returns the gradient of the cost function:\n",
    "\n",
    "$$\\nabla{C_x} = \\left[ \\frac{\\partial C_x}{\\partial b}, \\frac{\\partial C_x}{\\partial w} \\right]$$\n",
    "\n",
    "as a tuple $(\\nabla{b}, \\nabla{w})$. \n",
    "\n",
    "First we populate all $z$ and $a$ vectors for each layer during \"feedforward\" pass.\n",
    "\n",
    "Then we go backwards in the network starting at the last layer $L$ and calculate:\n",
    "\n",
    "$$\\delta_L = C^\\prime(z_L) $$\n",
    "\n",
    "This derivative of the cost function $C$ can be expressed using _chain rule_ as:\n",
    "\n",
    "$$ C^\\prime(a(z))a^\\prime(z) = \\frac{dC}{da(z)}a^\\prime(z)$$ \n",
    "\n",
    "note that\n",
    "\n",
    "$$ a(z) \\equiv \\sigma(z) $$\n",
    "\n",
    "and so $\\delta$ can be calculated as follows (this is the same as equation **BP1** in the book):\n",
    "\n",
    "$$ \\delta_L = C^\\prime(a_L) \\circ \\sigma^\\prime(z_L) $$\n",
    "\n",
    "where $\\circ$ denotes Hadamard product. After that we can calculate gradient vector elements for weights and biases: \n",
    "\n",
    "$$ \\nabla{b_L} = \\delta_{L} \\tag{equation BP3} $$\n",
    "\n",
    "$$ \\nabla{w_L} = \\delta_L \\cdot  a_{L-1}^\\top \\tag{equation BP4} $$\n",
    "\n",
    "Note that the original equation **BP4** was written as:\n",
    "\n",
    "$$ \\frac{\\partial{C}}{\\partial{w^L_{jk}}} = a^{L-1}_k\\delta^L_j $$\n",
    "\n",
    "but we want the shape of $\\nabla{w}$ to be the same as the shape of $w$, so we swap the operands while transposing the activation vector so it is now a row vector:\n",
    "\n",
    "$$ \\nabla{w} = \\begin{bmatrix} \n",
    "                \\delta_1 \\\\\n",
    "                \\delta_2 \\\\\n",
    "                \\vdots \\\\\n",
    "                \\delta_J \\\\\n",
    "               \\end{bmatrix} \\cdot \\begin{bmatrix} a_1 a_2 \\dots a_K \\end{bmatrix} $$\n",
    "\n",
    "where $J$ is the number of neurons in the $L$th layer, and $K$ is the number of neurons in the $(L-1)$th layer. This multiplication will result in the matrix of $(J, K)$ shape, that is with $J$ rows and $K$ columns. Remember that in the book the weight matrix elements $w^l_{jk}$ are defined as connection weights from $k$th neuron in the $l-1$ layer to the $j$th neuron in the $l$th layer, so the weight matrix would be the same $(J, K)$ shape.\n",
    "\n",
    "Then we recalculate $\\delta$ for each _preceeding_ layer $i$, starting with $i = L - 1$:\n",
    "\n",
    "$$ \\delta_{i} = w_{i+1}^\\top \\cdot \\delta_{i+1} \\circ \\sigma^\\prime(z_i) \\tag{equation BP2}$$\n",
    "\n",
    "and using $\\delta_i$ calculate new $\\nabla{b_i}$ and $\\nabla{w_i}$\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `nn`: the neural network object returned by the `init_network` function.\n",
    "- `mini_batch`: the mini-batch of tuples (x, y) where x is the input and y is the correct desired output\n",
    "\n",
    "Mini-batch is of the shape (($n$, 784, 1), ($n$, 10, 1)) where $n$ is the size (number of samples) in mini-batch.\n",
    "\n",
    "The function also calculates $z = w \\cdot x + b$ and the activation $a = \\sigma(z)$. $x, y, z, b, w$ are all `ndarrays`, and $w \\cdot x$ is the dot product.\n",
    "\n",
    "We're using `matmul` instead of `dot` to compute the dot product because we now have all matrices and vectors \"stacked\" in another dimension, each stacked \"layer\" corresponding to a sample in the mini-batch, and `matmul` does exactly what we want -- calculates dot products for each \"layer\" in the stack, and then stacks the results in the same manner.\n",
    "\n",
    "We're also using `transpose` method with axes $(0, 2, 1)$ to transpose only stacked 2D matrices, and not the whole 3D tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_backprop(nn, mini_batch):\n",
    "    nabla_b = [np.zeros(b.shape) for b in nn.biases]\n",
    "    nabla_w = [np.zeros(w.shape) for w in nn.weights]\n",
    "    \n",
    "    # mini-batch components\n",
    "    ax, ay = tuple(t for t in np.asarray(mini_batch).transpose())\n",
    "    \n",
    "    # ax and ay are one-dimensional (length of mini-batch) arrays of arrays of inputs x and outputs y\n",
    "    # we want to convert them to multi-dimensional arrays:\n",
    "    x = np.stack(ax)\n",
    "    y = np.stack(ay)  \n",
    "\n",
    "    # feedforward\n",
    "    activation = x    # first layer activation is just its input\n",
    "    activations = [x] # list to store all activations, layer by layer\n",
    "    zs = []           # list to store all z vectors, layer by layer\n",
    "\n",
    "    for b, w in zip(nn.biases, nn.weights):\n",
    "        z = np.matmul(w, activation) + b  # calculate z for the current layer\n",
    "        zs.append(z)                      # store\n",
    "        activation = sigmoid(z)           # layer output\n",
    "        activations.append(activation)    # store\n",
    "\n",
    "    # backward pass\n",
    "\n",
    "    # 1. starting from the output layer\n",
    "    delta = cost_derivative(activations[-1], y) * sigmoid_prime(zs[-1]) \n",
    "    nabla_b[-1] = delta.sum(axis = 0)\n",
    "    nabla_w[-1] = np.matmul(delta, activations[-2].transpose(0, 2, 1)).sum(axis = 0)\n",
    "    \n",
    "    # 2. continue back to the input layer (i is the layer index, we're using i instead of l\n",
    "    #    to improve readability -- l looks too much like 1)\n",
    "    for i in range(2, nn.num_layers): # starting from the next-to-last layer\n",
    "        z = zs[-i]\n",
    "        sp = sigmoid_prime(z)\n",
    "        delta = np.matmul(nn.weights[-i + 1].transpose(), delta) * sp\n",
    "        \n",
    "        nabla_b[-i] = delta.sum(axis = 0)\n",
    "        nabla_w[-i] = np.matmul(delta, activations[-i - 1].transpose(0, 2, 1)).sum(axis = 0)\n",
    "        \n",
    "    return (nabla_b, nabla_w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MNIST data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data load function\n",
    "\n",
    "Returns the MNIST data as a tuple containing the training data, the validation data, and the test data.\n",
    "\n",
    "The `training_data` is returned as a tuple with two entries. The first entry contains the actual training images.  This is an `ndarray` with 50,000 entries.  Each entry is, in turn, an `ndarray` with 784 values, representing the 28 * 28 = 784 pixels in a single MNIST image.\n",
    "\n",
    "The second entry in the `training_data` tuple is an `ndarray` containing 50,000 entries.  Those entries are just the digit values (0...9) for the corresponding images contained in the first entry of the tuple.\n",
    "\n",
    "The `validation_data` and `test_data` are similar, except each contains only 10,000 images.\n",
    "\n",
    "This is a nice data format, but for use in neural networks it's helpful to modify the format of the `training_data` a little.\n",
    "\n",
    "That's done in the wrapper function `load_data_wrapper()`, see below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data():\n",
    "    f = gzip.open('mnist.pkl.gz', 'rb')\n",
    "    training_data, validation_data, test_data = pickle.load(f, encoding=\"bytes\")\n",
    "    f.close()\n",
    "    \n",
    "    return (training_data, validation_data, test_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data reshaper function\n",
    "\n",
    "Returns a tuple containing `(training_data, validation_data, test_data)`. Based on `load_data`, but the format is more convenient for use in our implementation of neural networks.\n",
    "\n",
    "In particular, `training_data` is a list containing 50,000 2-tuples `(x, y)`.  `x` is a 784-dimensional `ndarray`\n",
    "containing the input image.  `y` is a 10-dimensional `ndarray` representing the unit vector corresponding to the\n",
    "correct digit for `x`.\n",
    "\n",
    "`validation_data` and `test_data` are lists containing 10,000 2-tuples `(x, y)` tuples. In each case, `x` is a 784-dimensional `ndarray` containing the input image, and `y` is the corresponding classification, i.e., the digit values (integers) corresponding to `x`.\n",
    "\n",
    "Obviously, this means we're using slightly different formats for the training data and the validation / test data.  These formats turn out to be the most convenient for use in our neural network code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_data_wrapper():\n",
    "    tr_d, va_d, te_d = load_data()\n",
    "    \n",
    "    training_inputs = [np.reshape(x, (784, 1)) for x in tr_d[0]]\n",
    "    training_results = [one_hot_encode(y) for y in tr_d[1]]\n",
    "    training_data = zip(training_inputs, training_results)\n",
    "    validation_inputs = [np.reshape(x, (784, 1)) for x in va_d[0]]\n",
    "    validation_data = zip(validation_inputs, va_d[1])\n",
    "    test_inputs = [np.reshape(x, (784, 1)) for x in te_d[0]]\n",
    "    test_data = zip(test_inputs, te_d[1])\n",
    "    \n",
    "    return (list(training_data), list(validation_data), list(test_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One-hot encoder\n",
    "\n",
    "Returns a 10-dimensional unit vector with a 1.0 in the $j$th position and zeroes elsewhere.  This is used to convert a digit (0...9) into a corresponding desired output from the neural network.\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `j`: the index in the array to set to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def one_hot_encode(j):\n",
    "    e = np.zeros((10, 1))\n",
    "    e[j] = 1.0\n",
    "    \n",
    "    return e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Utility function\n",
    "\n",
    "Prints the shape of a given `ndarray`. The first number is the number of rows, the second number is the number of columns.\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `name`: the name of the data\n",
    "- `data`: the `ndarray` of data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_shape(name, data):\n",
    "    print('Shape of {0}: {1}'.format(name, data.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main program\n",
    "\n",
    "Create and train a simple network to recognize handwritten digits based on MNIST data.\n",
    "\n",
    "We use 30 neurons in the intermediate layer and with the learning rate of 3.0, mini-batch size of 10, the network should achieve about 95% accuracy in 15 epochs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Network layer 2\n",
      "Shape of weights: (30, 784)\n",
      "Shape of biases: (30, 1)\n",
      "\n",
      "Network layer 3\n",
      "Shape of weights: (10, 30)\n",
      "Shape of biases: (10, 1)\n",
      "\n",
      "Learning process started...\n",
      "\n",
      "Epoch  1: accuracy 81.85%\n",
      "Epoch  2: accuracy 83.44%\n",
      "Epoch  3: accuracy 92.61%\n",
      "Epoch  4: accuracy 92.99%\n",
      "Epoch  5: accuracy 93.25%\n",
      "Epoch  6: accuracy 93.99%\n",
      "Epoch  7: accuracy 94.05%\n",
      "Epoch  8: accuracy 94.25%\n",
      "Epoch  9: accuracy 94.85%\n",
      "Epoch 10: accuracy 94.64%\n",
      "Epoch 11: accuracy 94.6%\n",
      "Epoch 12: accuracy 94.52%\n",
      "Epoch 13: accuracy 94.61%\n",
      "Epoch 14: accuracy 94.5%\n",
      "Epoch 15: accuracy 95.07%\n",
      "\n",
      "Learning process complete in 75 seconds (5.0 seconds per epoch)!\n",
      "\n",
      "Validation (with yet unseen data): accuracy 95.26%\n"
     ]
    }
   ],
   "source": [
    "training_data, validation_data, test_data = load_data_wrapper() # load data\n",
    "nn = init_network([784, 30, 10])\n",
    "\n",
    "for l in range(0, nn.num_layers - 1):\n",
    "    print('\\nNetwork layer {0}'.format(l + 2)) # disregard the input layer\n",
    "    print_shape('weights', nn.weights[l])\n",
    "    print_shape('biases', nn.biases[l])\n",
    "    \n",
    "# hyper parameters\n",
    "epochs = 15\n",
    "mini_batch_size = 10\n",
    "learning_rate = 3.0\n",
    "    \n",
    "print('\\nLearning process started...\\n')\n",
    "\n",
    "time_start = time()\n",
    "    \n",
    "learn(nn, training_data, epochs, mini_batch_size, learning_rate, test_data)\n",
    "\n",
    "time_end = time()\n",
    "\n",
    "time_elapsed = time_end - time_start\n",
    "\n",
    "print('\\nLearning process complete in {0} seconds ({1} seconds per epoch)!\\n'.format(f'{time_elapsed:.0f}', f'{time_elapsed / epochs:.1f}'))\n",
    "\n",
    "print('Validation (with yet unseen data): accuracy {0}%'.format(100.0 * evaluate(nn, validation_data) / len(validation_data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load image function\n",
    "\n",
    "This function loads an image from a file and returns an `ndarray` of floats with [784x1] dimensions, the same as in the MNIST set.\n",
    "\n",
    "The file is assumed to be 28x28 pixels, grayscale. _There is no error checking of any kind._\n",
    "\n",
    "Parameters:\n",
    "\n",
    "- `file_name`: the path to the image file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_image(file_name):\n",
    "    digit = Image.open(file_name)\n",
    "    \n",
    "    # invert, so that background is black (zeros) \n",
    "    digit = ImageOps.invert(digit)\n",
    "    \n",
    "    pixels = digit.load()\n",
    "    \n",
    "    return np.array(digit).reshape((784, 1)) / 255"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image recognition function\n",
    "\n",
    "This function takes a path template and the file name (which must be just a single-digit number corresponding to the digit image in the file), loads the image, and passes it through the neural network. The result is then compared with the file number.\n",
    "\n",
    "Parameters:\n",
    "- `path`: path template where the file name is inserted\n",
    "- `file`: the number corresponding to the file name e.g. 1 => '1.png'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def recognize_image(path, file):\n",
    "    x = load_image(path.format(file))\n",
    "\n",
    "    y = feedforward(nn, x)\n",
    "    \n",
    "    bitmap = x.reshape((28, 28))\n",
    "    \n",
    "    file_num = int(file)\n",
    "    result = y.argmax()\n",
    "    \n",
    "    if file_num == result:\n",
    "        ev = 'correctly'\n",
    "    else:\n",
    "        ev = 'incorrectly'\n",
    " \n",
    "    print(file, 'was', ev, 'recognized as', result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Real-life test!\n",
    "\n",
    "Load images of hand-written digits from our custom-made files and process them with our neural network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Non-MNIST digits:\n",
      "\n",
      "0 was incorrectly recognized as 3\n",
      "1 was correctly recognized as 1\n",
      "2 was correctly recognized as 2\n",
      "3 was correctly recognized as 3\n",
      "4 was correctly recognized as 4\n",
      "5 was correctly recognized as 5\n",
      "6 was correctly recognized as 6\n",
      "7 was correctly recognized as 7\n",
      "8 was correctly recognized as 8\n",
      "9 was correctly recognized as 9\n"
     ]
    }
   ],
   "source": [
    "print('Non-MNIST digits:\\n')\n",
    "\n",
    "for file in range(0,10):\n",
    "    recognize_image('./non-MNIST-digits/{0}.png', file)"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}