{"cells":[{"metadata":{"_cell_guid":"fc067bac-2132-47ad-aeb8-244e0faf0b1e","_uuid":"423a1ed04f75ef0903eb9661111d88c68bdbf6e7"},"cell_type":"markdown","source":"# Unsupervised learning with Self-Organizing Maps: MNIST\n\nIn this notebook, I will provide a short example of how Kohonen Self-Organizing Maps (SOM) can be used for dimensionality reduction and unsupervised learning. I will use the MNIST dataset provided by Kaggle to train a SOM and project the handwritten digits to a two-dimensional map that (hopefully) preserves the topological property of the original dataset. Although a number of alternative methods can be used to obtain competitive results in less time (in particular supervised learning with CNN is known to be extremely effective), the visual nature of the dataset and the easy interpretation of its elements (common handwritten digits) make it ideal to explore the capabilities of this method and understand how SOM work.\n\nThis notebook was inspired by [Interactive Intro to Dimensionality Reduction](https://www.kaggle.com/arthurtok/interactive-intro-to-dimensionality-reduction), by Anisotropic, where a number of common dimensionality reduction and clustering methods are presented. "},{"metadata":{"_cell_guid":"1ab74b89-2846-4902-aec9-44528d06d4b0","_uuid":"067cb3972bdb7cd0ba23a57a502e6b05ee14271a"},"cell_type":"markdown","source":"## What are SOM?\n\nIn brief, Self-organizing maps are a type of artificial neural network based on competitive learning (at variance to error-correcting learning typical of other NNs). The idea is to iteratively adapt a connected two-dimensional matrix of vectors (or nodes) to the higher-dimensional topology of the input dataset. \nAt each cycle, a node is selected and its elements (the weights) are updated, together with those of its neighbors, to approach a randomly chosen datapoint from the training set. The competitive element comes into play during the update stage, since the closest node (according to a chosen metric) to the extracted datapoint is selected for the weights update at each iteration.\n\nSOMs are particularly suited for cases where low-dimensional manifolds are hidden in higher dimensions and are often used together and/or competing with other dimensionality reduction methods and in particular Principal Component Analysis (PCA) for which it could be seen as a non-linear generalization: an exhaustive explanation of SOM's advantages and disadvantages, however, is beyond the scope of this notebook, but there are plenty of resources online for those who would like to know more."},{"metadata":{"_cell_guid":"e1ebe2c1-1113-4913-9017-4002bea5a3c0","_uuid":"8dbd88a992fd7ffbe16e5a22cb036b3abdc74224"},"cell_type":"markdown","source":"## Preparing the Data"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":false,"collapsed":true},"cell_type":"code","source":"#As usual we start importing a number of libraries that will be come in handy later on\nimport numpy as np \nimport pandas as pd \nimport seaborn as sns\nfrom imageio import imwrite\n#from scipy.misc import imsave\nimport matplotlib\nimport matplotlib.pyplot as plt\nfrom matplotlib.offsetbox import OffsetImage, AnnotationBbox\nfrom matplotlib.cbook import get_sample_data\n%matplotlib inline\n\nimport plotly.offline as py\npy.init_notebook_mode(connected=True)\nimport plotly.graph_objs as go\nimport plotly.tools as tls\n\nfrom sklearn.preprocessing import StandardScaler\nfrom PIL import Image, ImageChops","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"0aec18e5-3dea-4fa2-84a8-a58a545756e3","_uuid":"defe9ee1169f4686ac692ac777d6c3175f5b1976"},"cell_type":"markdown","source":"We will use here SimpSOM, a python library for SOM freely available on [PyPI](https://pypi.python.org/pypi/SimpSOM) and [Github](https://github.com/fcomitani/SimpSOM). "},{"metadata":{"collapsed":true,"_cell_guid":"2af08c97-9730-4a59-abc8-35cf9164ef5b","_uuid":"f172b6816ebdccab0332d56c7fab2d940856f064","trusted":false},"cell_type":"code","source":"import SimpSOM as sps","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"229fb11e-298c-4c64-91d5-6aae11af54c1","_uuid":"4e8c398d08f0e4d2ed2af5d1ab7c3437518b94c6"},"cell_type":"markdown","source":"To help the convergence of the map it is a good idea to limit the number of points on which the training will be done. We should thus choose a limited number of landmark points that possibly well represent the distribution of the entire population in the high dimensionality space. In our case we can assume that the distribution is uniform across the ten digits, we can then pick randomly a subset of these for the training. It is important to note, however, that although a higher number of training points does, in theory, increase the accuracy of the mapping by taking into account more variability in the images, it could also make the mapping process more complicated and hinder its convergence. Sometimes you just need a little less data!   \n\nHere, we will select only 500 points. Let's check that indeed they are equally distributed. It's a good idea to also fix the random seed for reproducibility and debugging purposes since the final results of the mapping will be dependent on its value.\n\nAs shown in other kernels, the images data is unfolded into 28\\*28 arrays, with each element containing a single black and white pixel color value. We will keep it this way, and normalize the arrays to improve the map learning."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false,"collapsed":true},"cell_type":"code","source":"np.random.seed(0)\n\n#We load here the data from the provided training set, we randomly select 500 landmark points and separate the labels. \ntrain = pd.read_csv('../input/train.csv')\ntrain = train.sample(n=500, random_state=0)\nlabels = train['label']\ntrain = train.drop(\"label\",axis=1)\n\n#Let's plot the distribution and see if the distribution is uniform\nsns.distplot(labels.values,bins=np.arange(-0.5,10.5,1))\n\n#Then we normalize the data, a crucial step to the correct functioning of the SOM algorithm\ntrainSt = StandardScaler().fit_transform(train.values)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"a70c1545-461a-4e74-8a38-ef6d1ab6875b","_uuid":"1d8a89cddebbb14056eb17f51c2f5f4acdaa4c9a"},"cell_type":"markdown","source":"The resulting distribution is not great, but it's what we can expect by taking only 500 points. It will have to do.\n\n## Training the Map\n\nNow we can proceed in setting up a map and train it with our landmark points.\n\nFirst, we need to build a network, it can be as big as one wants, but it should contain at least enough nodes to map all our landmark points separated, possibly more. Obviously, a bigger map allows to distinguish more subtle topological features, but comes at a cost, as the computation time increases with each added node.\n\nAgain, given the limited resources available on interactive kernels, we will build a 40x40 map and we will train it for only 5000 epochs (a good rule of thumb is to cycle for at least 10 times the number of training points) with a relatively high initial learning rate of 0.1. We activate periodic boundary conditions, to avoid artifacts at the borders and we chose PCA as weights initialization.\nThe nodes weights can be randomly chosen either from a uniform distribution within the minimum and maximum values of our training dataset (good thing we normalized it beforehand!) or from the space spanned by the two first PCA vectors. We chose here the latter as it may help the map converge faster.\n\nWith these parameters, it will still take around 6 hours on Kaggle's CPU."},{"metadata":{"collapsed":true,"_cell_guid":"07a1ccf4-29c5-4f07-97cb-2bbeac2ea4c0","_uuid":"50137fa5ecec01ff788f2e5311dff3dfd3fe1970","trusted":false},"cell_type":"code","source":"#We build a 40x40 network and initialise its weights with PCA \nnet = sps.somNet(40, 40, trainSt, PBC=True, PCI=True)\n\n#Now we can train it with 0.1 learning rate for 10000 epochs\nnet.train(0.1, 10000)\n\n#We print to screen the map of the weights differences between nodes, this will help us identify cluster centers \nnet.diff_graph(show=True,printout=True)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"24e5ddd1-4ea5-4f93-bad5-845f81709776","_uuid":"6f6bd7fdf9b096de249f4067504c6c8653858692"},"cell_type":"markdown","source":"Nice! We now have a trained map, but we need to understand what it means. In the plot above the color, the scale goes from dark blue to yellow for nodes whose neighbors weights are less or more distant (or simply different) respectively. This means that darker areas represent cluster basins where all neighbor nodes are similar, while yellow areas are high-difference border regions between basins. Just think of mountains and valleys! \n\n## Visualizing the Results\n\nWhat we need to do next is simply project our data (any MNIST data) on the map and see where each image gets mapped.In our case, the interpretation of the data is relatively easy, since all our pictures are labeled with the corresponding category. In pure unsupervised learning projects, where labels are not known, this can be difficult and a simple visual inspection may be problematic. A good approach can be to extract the weighted vectors from the basin center nodes and have a look at them as 28x28 images. Here we can do it for a few basins and see that indeed they are archetypes of specific digits: it works!"},{"metadata":{"collapsed":true,"_cell_guid":"f3bce337-749d-4945-9a04-011efba1069d","_uuid":"ede70c3adc691db815517f6c601cd99e21cd98f6","trusted":false},"cell_type":"code","source":"#Here we first define a few useful functions\ndef autocrop(fileName):\n    im = Image.open(fileName)\n    im=im.crop((0,100,2900,im.size[1]))\n    bg = Image.new(im.mode, im.size, im.getpixel((0,0)))\n    diff = ImageChops.difference(im, bg)\n    diff = ImageChops.add(diff, diff, 2.0, -100)\n    bbox = diff.getbbox()\n    if bbox:\n        return im.crop(bbox)\n\ndef posMap(x,y):\n     if y%2==0:\n        return [y, x*2/np.sqrt(3)*3/4]\n     else:\n        return [y+0.5, x*2/np.sqrt(3)*3/4]\n    \ndef posCount(x,y):\n     return y*40+x\n\ndef posCountR(x):\n     return [np.int(x%40),np.int(x/40)]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e9d9ac5d-7ccf-4a5f-900b-b584c158e6b6","_uuid":"50c7fcf138c6d5e43e038fc2ed7d8939dd68a13c","scrolled":false,"trusted":false,"collapsed":true},"cell_type":"code","source":"#Let's print a few trained nodes' weights and see how good they are\nlistNodes=[[20,0],[23,11],[1,6],[13,37],[7,33],[18,31]]\nlistCount=[posCount(20,0), posCount(23,11), posCount(1,6), posCount(13,37), posCount(7,33), posCount(18,31)]\n\ni=0\nfor node in net.nodeList:\n    if i in listCount:\n        print('Node\\'s position: {:d} {:d}'.format(posCountR(i)[1], posCountR(i)[0]) )\n        plt.imshow(np.asarray(node.weights).reshape(28,28))\n        plt.axis('off')\n        plt.show()\n    i+=1","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"8cc9d1df-a9d9-4e16-ab70-9f668b2d8457","_uuid":"7f4d286e8791668a3d27610d00d7974c11ce3973"},"cell_type":"markdown","source":"As expected, some of them are quite similar to specific digits archetypes, while others still contain generic and undefined shapes, either because they lay closer to a border region between different categorical clusters, or simply because in need of more training.\nLet's now map our training points to see where they ended up by matching them with their closest node. We map the results interactively with plotly, hovering on each point with the mouse we can see the original label. To help visualize the clusters, each label has been colored differently in the scatter plot."},{"metadata":{"_cell_guid":"8519bd95-73fa-43d6-8183-c2e78c577404","_uuid":"8f6082d40e327945806abf2dc230809332a1dedb","trusted":false,"collapsed":true},"cell_type":"code","source":"projData=net.project(trainSt[:500])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2b799950-0071-425f-8c6d-f3ac217d38de","_uuid":"3b26f212682d9100bdf4aec0213def12f006f8ad","trusted":false,"collapsed":true},"cell_type":"code","source":"#We first save a cropped version of the original map to superimpose and then we add the scatterpoints\ncropped = autocrop('nodesDifference.png')\ncropped.save('cropped.png')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"6a8e355e-37e9-4883-97df-281d066a14a8","_uuid":"299c68872b8aab2b484aafe7469c87586d0af234","trusted":false,"collapsed":true},"cell_type":"code","source":"#And here we prepare the plotly graph. \ntrace0 = go.Scatter(\n    x = [x for x,y in projData],\n    y = [y for x,y in projData],\n#    name = labels,\n    hoveron = [str(n) for n in labels],\n    text = [str(n) for n in labels],\n    mode = 'markers',\n    marker = dict(\n        size = 8,\n        color = labels,\n        colorscale ='Jet',\n        showscale = False,\n        opacity = 1\n    ),\n    showlegend = False\n\n)\ndata = [trace0]\n\nlayout = go.Layout(\n    images= [dict(\n                  source= \"cropped.png\",\n                  xref= \"x\",\n                  yref= \"y\",\n                  x= -0.5,\n                  y= 39.5*2/np.sqrt(3)*3/4,\n                  sizex= 40.5,\n                  sizey= 40*2/np.sqrt(3)*3/4,\n                  sizing= \"stretch\",\n                  opacity= 0.5,\n                  layer= \"below\")],\n    width = 800,\n    height = 800,\n    hovermode= 'closest',\n    xaxis= dict(\n        range=[-1,41],\n        zeroline=False,\n        showgrid=False,\n        ticks='',\n        showticklabels=False\n    ),\n    yaxis=dict(\n        range=[-1,41],\n        zeroline=False,\n        showgrid=False,\n        ticks='',\n        showticklabels=False\n    ),\n    showlegend= True\n)\n\n\nfig = dict(data=data, layout=layout)\npy.iplot(fig, filename='styled-scatter')","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"85cf53d7-dea2-4d26-bbf3-ec15731d75a8","_uuid":"3054737e8bfea4b8c1a897c2d81614f60451738e"},"cell_type":"markdown","source":"Nice! The algorithm is able to clearly separate in clusters some of our categories. Their centers are mostly found in the \"basins\" where the weights difference is minimal (in dark blue). The topology of the high dimensionality space is also somehow respected, with similar shapes mapped closed together (e.g. handwritten 7 and 9) and others far apart (such 1 and 0).  Unfortunately, however, the accuracy is not great, as there are a number of outliers. This issue could be easily solved by fine-tuning the algorithm parameters, with further training and a better representative set of landmark points. Introducing more variability in the training set, by adding random rotations and translations, could really make the difference."},{"metadata":{"_cell_guid":"8ad057ff-1c63-48a1-a21a-3410f8b84b84","_uuid":"907ea646a0ca834c03e56e8c9da176d37c94d901"},"cell_type":"markdown","source":"\n\n## Conclusion\n\n"},{"metadata":{"collapsed":true,"_cell_guid":"a7393c98-57ce-4d3f-acb5-185a77c1c3ac","_uuid":"8da6d18f324f9e18353ea7d264d92e9b302ff912"},"cell_type":"markdown","source":"Self-Organizing Maps represent an interesting alternative to more commonly used unsupervised learning and dimensionality reduction algorithms. Visualizing the topology of the dataset can be helpful in identifying hidden patterns and relationships between clusters of datapoints. The algorithm can be easily extended and used in conjunction with other methods (such as PCA to initialize the weights or k-means to cluster the projected data on top of the map) to tackle more challenging problems. "}],"metadata":{"language_info":{"version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"version":3,"name":"ipython"},"pygments_lexer":"ipython3","nbconvert_exporter":"python","name":"python","file_extension":".py"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"nbformat":4,"nbformat_minor":1}