Initial files

.gitignore

@@ -0,0 +1 @@
+__pycache__

expand_mnist.py

@@ -0,0 +1,60 @@
+"""expand_mnist.py
+~~~~~~~~~~~~~~~~~~
+
+Take the 50,000 MNIST training images, and create an expanded set of
+250,000 images, by displacing each training image up, down, left and
+right, by one pixel.  Save the resulting file to
+../data/mnist_expanded.pkl.gz.
+
+Note that this program is memory intensive, and may not run on small
+systems.
+
+"""
+
+from __future__ import print_function
+
+#### Libraries
+
+# Standard library
+import cPickle
+import gzip
+import os.path
+import random
+
+# Third-party libraries
+import numpy as np
+
+print("Expanding the MNIST training set")
+
+if os.path.exists("../data/mnist_expanded.pkl.gz"):
+    print("The expanded training set already exists.  Exiting.")
+else:
+    f = gzip.open("../data/mnist.pkl.gz", 'rb')
+    training_data, validation_data, test_data = cPickle.load(f)
+    f.close()
+    expanded_training_pairs = []
+    j = 0 # counter
+    for x, y in zip(training_data[0], training_data[1]):
+        expanded_training_pairs.append((x, y))
+        image = np.reshape(x, (-1, 28))
+        j += 1
+        if j % 1000 == 0: print("Expanding image number", j)
+        # iterate over data telling us the details of how to
+        # do the displacement
+        for d, axis, index_position, index in [
+                (1,  0, "first", 0),
+                (-1, 0, "first", 27),
+                (1,  1, "last",  0),
+                (-1, 1, "last",  27)]:
+            new_img = np.roll(image, d, axis)
+            if index_position == "first": 
+                new_img[index, :] = np.zeros(28)
+            else: 
+                new_img[:, index] = np.zeros(28)
+            expanded_training_pairs.append((np.reshape(new_img, 784), y))
+    random.shuffle(expanded_training_pairs)
+    expanded_training_data = [list(d) for d in zip(*expanded_training_pairs)]
+    print("Saving expanded data. This may take a few minutes.")
+    f = gzip.open("../data/mnist_expanded.pkl.gz", "w")
+    cPickle.dump((expanded_training_data, validation_data, test_data), f)
+    f.close()

mnist_average_darkness.py

@@ -0,0 +1,64 @@
+"""
+mnist_average_darkness
+~~~~~~~~~~~~~~~~~~~~~~
+
+A naive classifier for recognizing handwritten digits from the MNIST
+data set.  The program classifies digits based on how dark they are
+--- the idea is that digits like "1" tend to be less dark than digits
+like "8", simply because the latter has a more complex shape.  When
+shown an image the classifier returns whichever digit in the training
+data had the closest average darkness.
+
+The program works in two steps: first it trains the classifier, and
+then it applies the classifier to the MNIST test data to see how many
+digits are correctly classified.
+
+Needless to say, this isn't a very good way of recognizing handwritten
+digits!  Still, it's useful to show what sort of performance we get
+from naive ideas."""
+
+#### Libraries
+# Standard library
+from collections import defaultdict
+
+# My libraries
+import mnist_loader
+
+def main():
+    training_data, validation_data, test_data = mnist_loader.load_data()
+    # training phase: compute the average darknesses for each digit,
+    # based on the training data
+    avgs = avg_darknesses(training_data)
+    # testing phase: see how many of the test images are classified
+    # correctly
+    num_correct = sum(int(guess_digit(image, avgs) == digit)
+                      for image, digit in zip(test_data[0], test_data[1]))
+    print "Baseline classifier using average darkness of image."
+    print "%s of %s values correct." % (num_correct, len(test_data[1]))
+
+def avg_darknesses(training_data):
+    """ Return a defaultdict whose keys are the digits 0 through 9.
+    For each digit we compute a value which is the average darkness of
+    training images containing that digit.  The darkness for any
+    particular image is just the sum of the darknesses for each pixel."""
+    digit_counts = defaultdict(int)
+    darknesses = defaultdict(float)
+    for image, digit in zip(training_data[0], training_data[1]):
+        digit_counts[digit] += 1
+        darknesses[digit] += sum(image)
+    avgs = defaultdict(float)
+    for digit, n in digit_counts.iteritems():
+        avgs[digit] = darknesses[digit] / n
+    return avgs
+
+def guess_digit(image, avgs):
+    """Return the digit whose average darkness in the training data is
+    closest to the darkness of ``image``.  Note that ``avgs`` is
+    assumed to be a defaultdict whose keys are 0...9, and whose values
+    are the corresponding average darknesses across the training data."""
+    darkness = sum(image)
+    distances = {k: abs(v-darkness) for k, v in avgs.iteritems()}
+    return min(distances, key=distances.get)
+
+if __name__ == "__main__":
+    main()

mnist_loader.py

@@ -0,0 +1,78 @@
+# %load mnist_loader.py
+"""
+mnist_loader
+~~~~~~~~~~~~
+A library to load the MNIST image data.  For details of the data
+structures that are returned, see the doc strings for ``load_data``
+and ``load_data_wrapper``.  In practice, ``load_data_wrapper`` is the
+function usually called by our neural network code.
+"""
+
+#### Libraries
+# Standard library
+import pickle
+import gzip
+
+# Third-party libraries
+import numpy as np
+
+def load_data():
+    """Return the MNIST data as a tuple containing the training data,
+    the validation data, and the test data.
+    The ``training_data`` is returned as a tuple with two entries.
+    The first entry contains the actual training images.  This is a
+    numpy ndarray with 50,000 entries.  Each entry is, in turn, a
+    numpy ndarray with 784 values, representing the 28 * 28 = 784
+    pixels in a single MNIST image.
+    The second entry in the ``training_data`` tuple is a numpy 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.
+    The ``validation_data`` and ``test_data`` are similar, except
+    each contains only 10,000 images.
+    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.
+    That's done in the wrapper function ``load_data_wrapper()``, see
+    below.
+    """
+    f = gzip.open('mnist.pkl.gz', 'rb')
+    training_data, validation_data, test_data = pickle.load(f, encoding="latin1")
+    f.close()
+    return (training_data, validation_data, test_data)
+
+def load_data_wrapper():
+    """Return 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.
+    In particular, ``training_data`` is a list containing 50,000
+    2-tuples ``(x, y)``.  ``x`` is a 784-dimensional numpy.ndarray
+    containing the input image.  ``y`` is a 10-dimensional
+    numpy.ndarray representing the unit vector corresponding to the
+    correct digit for ``x``.
+    ``validation_data`` and ``test_data`` are lists containing 10,000
+    2-tuples ``(x, y)``.  In each case, ``x`` is a 784-dimensional
+    numpy.ndarry containing the input image, and ``y`` is the
+    corresponding classification, i.e., the digit values (integers)
+    corresponding to ``x``.
+    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."""
+    tr_d, va_d, te_d = load_data()
+    training_inputs = [np.reshape(x, (784, 1)) for x in tr_d[0]]
+    training_results = [vectorized_result(y) for y in tr_d[1]]
+    training_data = zip(training_inputs, training_results)
+    validation_inputs = [np.reshape(x, (784, 1)) for x in va_d[0]]
+    validation_data = zip(validation_inputs, va_d[1])
+    test_inputs = [np.reshape(x, (784, 1)) for x in te_d[0]]
+    test_data = zip(test_inputs, te_d[1])
+    return (training_data, validation_data, test_data)
+
+def vectorized_result(j):
+    """Return a 10-dimensional unit vector with a 1.0 in the jth
+    position and zeroes elsewhere.  This is used to convert a digit
+    (0...9) into a corresponding desired output from the neural
+    network."""
+    e = np.zeros((10, 1))
+    e[j] = 1.0
+    return e

mnist_svm.py

@@ -0,0 +1,28 @@
+"""
+mnist_svm
+~~~~~~~~~
+
+A classifier program for recognizing handwritten digits from the MNIST
+data set, using an SVM classifier."""
+
+#### Libraries
+# My libraries
+import mnist_loader 
+
+# Third-party libraries
+from sklearn import svm
+
+def svm_baseline():
+    training_data, validation_data, test_data = mnist_loader.load_data()
+    # train
+    clf = svm.SVC()
+    clf.fit(training_data[0], training_data[1])
+    # test
+    predictions = [int(a) for a in clf.predict(test_data[0])]
+    num_correct = sum(int(a == y) for a, y in zip(predictions, test_data[1]))
+    print "Baseline classifier using an SVM."
+    print "%s of %s values correct." % (num_correct, len(test_data[1]))
+
+if __name__ == "__main__":
+    svm_baseline()
+

network.py

@@ -0,0 +1,149 @@
+# %load network.py
+
+"""
+network.py
+~~~~~~~~~~
+IT WORKS
+
+A module to implement the stochastic gradient descent learning
+algorithm for a feedforward neural network.  Gradients are calculated
+using backpropagation.  Note that I have focused on making the code
+simple, easily readable, and easily modifiable.  It is not optimized,
+and omits many desirable features.
+"""
+
+#### Libraries
+# Standard library
+import random
+
+# Third-party libraries
+import numpy as np
+
+class Network(object):
+
+    def __init__(self, sizes):
+        """The list ``sizes`` contains the number of neurons in the
+        respective layers of the network.  For example, if the list
+        was [2, 3, 1] then it would be a three-layer network, with the
+        first layer containing 2 neurons, the second layer 3 neurons,
+        and the third layer 1 neuron.  The biases and weights for the
+        network are initialized randomly, using a Gaussian
+        distribution with mean 0, and variance 1.  Note that the first
+        layer is assumed to be an input layer, and by convention we
+        won't set any biases for those neurons, since biases are only
+        ever used in computing the outputs from later layers."""
+        self.num_layers = len(sizes)
+        self.sizes = sizes
+        self.biases = [np.random.randn(y, 1) for y in sizes[1:]]
+        self.weights = [np.random.randn(y, x)
+                        for x, y in zip(sizes[:-1], sizes[1:])]
+
+    def feedforward(self, a):
+        """Return the output of the network if ``a`` is input."""
+        for b, w in zip(self.biases, self.weights):
+            a = sigmoid(np.dot(w, a)+b)
+        return a
+
+    def SGD(self, training_data, epochs, mini_batch_size, eta,
+            test_data=None):
+        """Train the neural network using mini-batch stochastic
+        gradient descent.  The ``training_data`` is a list of tuples
+        ``(x, y)`` representing the training inputs and the desired
+        outputs.  The other non-optional parameters are
+        self-explanatory.  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 substantially."""
+
+        training_data = list(training_data)
+        n = len(training_data)
+
+        if test_data:
+            test_data = list(test_data)
+            n_test = len(test_data)
+
+        for j in range(epochs):
+            random.shuffle(training_data)
+            mini_batches = [
+                training_data[k:k+mini_batch_size]
+                for k in range(0, n, mini_batch_size)]
+            for mini_batch in mini_batches:
+                self.update_mini_batch(mini_batch, eta)
+            if test_data:
+                print("Epoch {} : {} / {}".format(j,self.evaluate(test_data),n_test));
+            else:
+                print("Epoch {} complete".format(j))
+
+    def update_mini_batch(self, mini_batch, eta):
+        """Update the network's weights and biases by applying
+        gradient descent using backpropagation to a single mini batch.
+        The ``mini_batch`` is a list of tuples ``(x, y)``, and ``eta``
+        is the learning rate."""
+        nabla_b = [np.zeros(b.shape) for b in self.biases]
+        nabla_w = [np.zeros(w.shape) for w in self.weights]
+        for x, y in mini_batch:
+            delta_nabla_b, delta_nabla_w = self.backprop(x, y)
+            nabla_b = [nb+dnb for nb, dnb in zip(nabla_b, delta_nabla_b)]
+            nabla_w = [nw+dnw for nw, dnw in zip(nabla_w, delta_nabla_w)]
+        self.weights = [w-(eta/len(mini_batch))*nw
+                        for w, nw in zip(self.weights, nabla_w)]
+        self.biases = [b-(eta/len(mini_batch))*nb
+                       for b, nb in zip(self.biases, nabla_b)]
+
+    def backprop(self, x, y):
+        """Return a tuple ``(nabla_b, nabla_w)`` representing the
+        gradient for the cost function C_x.  ``nabla_b`` and
+        ``nabla_w`` are layer-by-layer lists of numpy arrays, similar
+        to ``self.biases`` and ``self.weights``."""
+        nabla_b = [np.zeros(b.shape) for b in self.biases]
+        nabla_w = [np.zeros(w.shape) for w in self.weights]
+        # feedforward
+        activation = x
+        activations = [x] # list to store all the activations, layer by layer
+        zs = [] # list to store all the z vectors, layer by layer
+        for b, w in zip(self.biases, self.weights):
+            z = np.dot(w, activation)+b
+            zs.append(z)
+            activation = sigmoid(z)
+            activations.append(activation)
+        # backward pass
+        delta = self.cost_derivative(activations[-1], y) * \
+            sigmoid_prime(zs[-1])
+        nabla_b[-1] = delta
+        nabla_w[-1] = np.dot(delta, activations[-2].transpose())
+        # Note that the variable l in the loop below is used a little
+        # differently to the notation in Chapter 2 of the book.  Here,
+        # l = 1 means the last layer of neurons, l = 2 is the
+        # second-last layer, and so on.  It's a renumbering of the
+        # scheme in the book, used here to take advantage of the fact
+        # that Python can use negative indices in lists.
+        for l in range(2, self.num_layers):
+            z = zs[-l]
+            sp = sigmoid_prime(z)
+            delta = np.dot(self.weights[-l+1].transpose(), delta) * sp
+            nabla_b[-l] = delta
+            nabla_w[-l] = np.dot(delta, activations[-l-1].transpose())
+        return (nabla_b, nabla_w)
+
+    def evaluate(self, test_data):
+        """Return 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."""
+        test_results = [(np.argmax(self.feedforward(x)), y)
+                        for (x, y) in test_data]
+        return sum(int(x == y) for (x, y) in test_results)
+
+    def cost_derivative(self, output_activations, y):
+        """Return the vector of partial derivatives \partial C_x /
+        \partial a for the output activations."""
+        return (output_activations-y)
+
+#### Miscellaneous functions
+def sigmoid(z):
+    """The sigmoid function."""
+    return 1.0/(1.0+np.exp(-z))
+
+def sigmoid_prime(z):
+    """Derivative of the sigmoid function."""
+    return sigmoid(z)*(1-sigmoid(z))