Kohonen Map with Numba

Kohonen.ipynb

{
 "metadata": {
  "name": "Kohonen Map"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# A Kohonen Map in Python (optimized by Numba)  \n",
      "\n",
      "Below is an implementation of a [Kohonen map](http://en.wikipedia.org/wiki/Self-organizing_map) (also known as a self-organizing map).   \n",
      "Briefly, it's a technique for mapping a high-dimensional space onto a much lower one. You might say, \"Well, principal components analysis can do that! Why do I need anything else?\"  \n",
      "PCA does find a projection from a high-D space down to a simpler one, but it does so linearly. Meaning, it can only rotate, scale, add and subtract dimensions. Self-organizing maps are free to do much weirder transformations when compressiong dimensionality. \n",
      "\n",
      "I needed a Kohonen map for my research, so I wrote this, and now I'm sharing it with you.\n",
      "\n",
      "Enjoy!  \n",
      "- Alex Wiltschko"
     ]
    },

    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from numba.decorators import jit as jit\n",
      "np.random.seed(1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 197
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Define helper functions\n",
      "\n",
      "These guys below are the core of the Kohonen map algorithm. They've been numba-ized for speed.  \n",
      "What's interesting is that they're still relatively readable (quite readable for how fast they are)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def scale(start, end, position_between_zero_and_one):\n",
      "    return double(start + position_between_zero_and_one * (end - start))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 198
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "nd0 = numba.double\n",
      "nd1 = numba.double[:]\n",
      "nd2 = numba.double[:,:]\n",
      "\n",
      "@jit(arg_types=[ nd2, nd1, nd1 ])\n",
      "def find_winner(som_weights, this_sample, abs_weight_distances):\n",
      "    winner_idx = 0\n",
      "    winner_distance = 1.0e100\n",
      "    num_nodes, num_dims = som_weights.shape\n",
      "\n",
      "    for i in range(num_nodes):\n",
      "        abs_weight_distances[i] = 0.0\n",
      "        for j in range(num_dims):\n",
      "            abs_weight_distances[i] += (((this_sample[j] - som_weights[i,j]) ** 2.0) ** 0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 208
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "@jit(arg_types=[ nd1, nd2, nd0, nd2, nd0, nd0 ])\n",
      "def update_weights(this_sample, som_weights, winner_idx, grid_indices, learning_rate, learning_spread):\n",
      "    \n",
      "    winner_x = grid_indices[winner_idx,0]\n",
      "    winner_y = grid_indices[winner_idx,1]\n",
      "    \n",
      "    num_nodes, num_dims = som_weights.shape\n",
      "    for i in range(num_nodes):\n",
      "        grid_distance = ((winner_x - grid_indices[i,0]) ** 2.0) + ((winner_y - grid_indices[i,1])**2.0)\n",
      "        dampening = e ** (-1.0 * grid_distance / ( 2.0 * learning_spread**2.0))\n",
      "        dampening *= learning_rate\n",
      "        \n",
      "        for j in range(num_dims):\n",
      "            som_weights[i,j] += dampening * (this_sample[j] - som_weights[i,j])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": []
      }
     ],
     "prompt_number": 207
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Setup parameters and allocate arrays"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# Generate some random colors for good times sake and what not have a good time yeah?\n",
      "X = np.random.random((16,3))\n",
      "\n",
      "# Some initial logistics\n",
      "n_samples, n_dims = X.shape\n",
      "\n",
      "# Construct the grid indices (2D for now)\n",
      "# The grid indices are stored in a (numNodePoints x numDimensionsToGrid) array. \n",
      "# It's easier to handle a list of the grid indices than it is to do anything otherwise. \n",
      "grid_size = (20,20)\n",
      "num_nodes = grid_size[0]*grid_size[1]\n",
      "raw_grid = np.mgrid[0:grid_size[0], 0:grid_size[1]]\n",
      "grid_indices = np.zeros((num_nodes, len(grid_size)), dtype='d')\n",
      "grid_indices[:,0] = raw_grid[0].ravel()\n",
      "grid_indices[:,1] = raw_grid[1].ravel()\n",
      "\n",
      "# Set parameters\n",
      "num_iterations = 200\n",
      "learning_rate_initial = 0.5\n",
      "learning_rate_final = 0.1\n",
      "learning_spread_initial = np.max(grid_size) / 5.0\n",
      "learning_spread_final = 1.0\n",
      "\n",
      "# Allocate the weight distances\n",
      "abs_weight_distances = np.zeros((num_nodes,), dtype='d')\n",
      "\n",
      "# Initialize the som_weights\n",
      "som_weights = np.random.random((num_nodes, n_dims))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 202
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Show off the initialized map"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "the_som = np.reshape(som_weights, (grid_size[0], grid_size[1], n_dims))\n",
      "figure(figsize=(5,5))\n",
      "imshow(the_som);\n",
      "\n",
      "figure(figsize=(4,4))\n",
      "axis('off')\n",
      "numRows = 4\n",
      "imshow(X[:,np.newaxis,:].reshape(numRows,X.shape[0]/numRows,3))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 203,
       "text": [
        "<matplotlib.image.AxesImage at 0x11630db10>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEzCAYAAAC7cS8aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtclVW+BvAHEBVRBAUqEFMEE9RARM07Yqko2dVOeaZM\nizqO03G6OdWMk+WkqcdKs5zJmXK0punUdJkzTYYxeL8leAVNQCjAC15QVMALcv7Ybo/xrg3sX7rX\nOi/P9x8/ay9+rCVuF+t99rvX9qqtra0FEZGNeOueABHR1caFjYhshwsbEdkOFzYish0ubERkO1zY\niMh2ml3rATIyMq71EETURA0fPlz5+DVf2ABgftoot2seDR4qGmv73jtEdX+d00lUd2aX7DbAO5bv\nEdUl/eo/RHVHH+ktqnvtYLKorsWNp0R1b2XNEdUVNf+VqO6D6BBR3d/X9xLV3ZSyW1Q3bdoWUd09\nsetEdV+tniWqu7XlRFFdxkfN3a7puGKbyz6tl6LHqnhvsEppzX7dUzDS9h0bdU/BOHkHTuqegpHE\nO7bc3FwsWbIENTU1SElJQUpKytWcFxGRmHjHtnTpUqSlpWH69On4+uuvUVFR4fb3aO/nJR3e1sJ9\nInVPwUjxcf11T8E40WFtdU/BSKKFrbKyEgAQGxuLkJAQxMXFIT8//6pOjIhISrSw5efnIzw8/HI7\nPDwceXl5bn8fZmxqzNjUmLFZMWNT88iroseqai9fdjoXs4baTrsrygEAPQKCGtUuuuBYYDs1i3ar\nDXQCAFTu2wUAaNW1Z6PaZ0sd7Rbh7rWdP3rnIua8/Cyt2Y+jNQd+1L6yP6fQ8SpX986D3Gpfd2nU\nfRuqAABdB/g1ql219QAAwC8xzK12ixvbAADOrC8DAPgPDG1U27l4OS87G9sO7OP4++3d7vh+3eJD\nG9Uu31gCAAjq38GtNuB4VXRNzl4AwJDu3RrVrl5bDABoOTjCrbaTcyFzXoKWHD39o3bd/g2Xnt4D\nHE/3Rred9h5xPP+6hUQ2qr3x8AYAQP/rBrjVBpIc36/W8Xzt5jWoUW1XvCTHFlVWVmLGjBmYO3cu\nAODdd99FfHw8EhISLF+bkZHB2z0UeLuHGm/3UOPtHlYdV2xzeR+b6FK0VatWAByvjJaVlWHnzp2I\nioqSfCsioqtOfCn68MMPY8mSJbhw4QJSUlIQEBDg9ve48hKV/k9pzX6+MqqwfcdGvjJaR96Bk3xl\nVEG8sMXGxuL111+/mnMhIroqtL7zgLs1Ne7W1Lhbs+JuTc0jr4o+XnCT2zVPlf+naKzju2T/0Hln\nJojqwn55SFSXk/OyqK7L+ThR3crP9orqfjE1R1T3yea3RXURvp+I6h65c0PDX6RwuHV4w1+kkPm0\n7O93ptl8Ud3SvWNEde/2GCCq+0uZ7MWDjLeTRHWP/XOc2zUr4HosrTu23avO6BzeWCc2F+iegpE2\nb+dN4HWdOrFd9xSMxPPYiMh2tC5sPZL8dQ5vrMB+XXRPwUj94nlLUV1tAuN1T8FI3LERke0wYzMQ\nMzY1ZmxWzNjUuGMjItthxmYgZmxqzNismLGpccdGRLbDjM1AzNjUmLFZMWNT446NiGyHGZuBmLGp\nMWOzYsamxh0bEdkOMzYDMWNTY8ZmxYxNzSOne5S9Ga18vDzvKMp2BSv79nl/Jhrrje6jRXX3/7vs\nKPIlcbKjkH/5uOuDOS+UB6JZUKKy79TEC6Lx0j+9QVS3fZXsCPOH+n4nqvv8+Bsu+/ZWrkV178HK\nvin7QkXj5S6sEtXteegLUd3dtbKj1r/54V7l49VVRThX4foY9tIHZUef9/i97N+943nZpfGRmY+6\nXZMF17/otO7YbopWL2pNnatFranrNli9qDVlLf066Z6CkZixEZHtaF3Yvss7qnN4Y10o36p7Ckba\nu3at7ikYp7qqSPcUjMQdGxHZDjM2AzFjU2PGZsWMTY07NiKyHWZsBmLGpsaMzYoZmxp3bERkO8zY\nDMSMTY0ZmxUzNjXu2IjIdpixGYgZmxozNitmbGrcsRGR7TBjMxAzNjVmbFbM2NQ8crrH4HHq0z3q\n83rpLtFYFVkPiOrufWiLqO7zZzeI6jqW+onq8pd3E9XVZhwW1X3w5Z9Fdde/niaq+/j9LFHdwu5t\nRXV37Gkpqnv72d+K6rzLnxLV3Z7r+jSY+qx/5hVR3R/+sVBUV/DoeFFdWfLVjaW07ti2bPhB5/DG\nOp0tO/LH7o7v49ljdZWsydE9BSMxYyMi29G6sPUd0FHn8MZqnXCT7ikYqV1Xnu9fV4ch3XVPwUjc\nsRGR7TBjMxAzNjVmbFbM2NS4YyMi22HGZiBmbGrM2KyYsalxx0ZEtsOMzUDM2NSYsVkxY1Pjjo2I\nbIcZm4GYsakxY7NixqbGHRsR2Q4zNgMxY1NjxmbFjE3NI6d7pO9Wn7yQV9gaJwLUff333Cwa680Z\n14vqFh+eKKpbWfqSqG7NxE0u+/LOlCDaX/07Z2NgJ9F4aYePi+oqohaL6j75k+ypVdWpj8u+s21K\nUdUvT9m3dtLnovEODGktqusRNENUN6uP7NSMh7erTz2pzL+I6wJc7086PbVTNN6/Lsh+LilvyX4p\nZ+3b5nZNKDq77NO6Y4uOcz2xpizan9mjSmC/cN1TME7nBKZJKvypEJHtaF3Y8nYU6hzeWHlnmD2q\nnNhcqnsKxinMvqh7Ckbijo2IbIcZm4GYsakxY7NixqbGnwoR2Q4zNgMxY1NjxmbFjE2NOzYish1m\nbAZixqbGjM2KGZsafypEZDvM2AzEjE2NGZsVMzY17tiIyHaYsRmIGZsaMzYrZmxqHjndY8DYJLdr\nygdMF401euRcUV1YVi9R3fQHD4vqhuw+Jar74M0SUV30kzeK6pov8xLV7Y/4SFT35LhVorqSrTeI\n6r4Z8V+iuj5b14jq2mW0F9XN3Co7pWN8t49FdegzRlTWd9WDorrDr7j/fPn2d8+57NO63GdddP+o\nkqagYM0e3VMw0u7TskXdzi7u26J7CkbiPpaIbEd8KTplyhT4+fnB29sbPj4+mD17ttvfo7e37PLP\n7roMidE9BSP1aN1B9xSM4921r+4pGOknZWwzZsxA69aykzaJiK6Vn3QpWltb+5MGZ8amxoxNjRmb\nFTM2NfGOzcvLCy+//DJCQ0MxbNgwJCYmXs15ERGJiRe2mTNnIigoCCUlJZgzZw6ioqIQGBjo1vdg\nxqbGjE2NGZsVMzY18cIWFBQEAOjQoQMSExOxdetW3Hrrrcqvzbq47fIi5rz8bKgdeal2R3k5ACDu\n0ngNtfeUOT7NJya0t1ttwDF+cfZqAEBEwtBGtY9VbQAAtPcb4FbbyXnZ6VzMGmqvKtgPAEjqEulW\n26nqsOMtOH7XeTeqXbnd0W4V714bEY4/Dm3eDQC4vl+PRrWdl5vORayxbeev1LzVjvsDo4e2aVT7\nRLnjI/0Cg+LdajsVX3C8LTCiWedGtSuyjgEAAnq3d6vt5Lz0dC5oDbUPbneMf0N8Z7fa3S59WFjB\n2q0AgC6DExvVvrjL8XZA754d3Wo7ra1w/H8eHBDUqLYrXrWCoOzs2bO4ePEi/Pz8UFFRgRkzZuCF\nF15AcHCw5WszMjLQNrWl8vtcueDVJb1B92Co7Abd7VNlu8c88Q26GS77Ctbscblre1p6g+62GaK6\n5o/KbtDtI7xB9556btDdfbrE5a6tZOta0Xgfe/gG3TUZskxsx1b1PC/u21Lvrk16g2434Q26L66S\n/bsf/l222zXf/u45DB8+XNkn2rGdPHkS8+bNAwC0adMGqampykWNiEgH0cIWGhp6eWH7KZixqTFj\nU2PGZsWMTY3vPCAi2+F7RQ3E+9jUeB+bFe9jU/PI6R4paf+mfPx8yVn4dmih7Js/4W3RWJHzZO+E\n2Oq/XFS3tzxNVLf+N9e77DtffBa+6eqfy/TUUaLxvjr0C1HdyvTxorpj/zlYVBdyYpHLvrbVvgi5\n0E3Zd/f8N0TjnXjfX1RXNfRDUd3KXf8S1Z3puFv5+OoTZRjasdhl3ZEPXb9IVZ/bJtwuqouIlj1f\nHlqT1fAX1fFEPa/faN2xuVrUmjrfCP5cVBJaqhe1pmzozXG6p2AkZmxEZDtaF7bzJWd1Dm+s88X8\nuahkV+/VPQXjrN65Q/cUjMQdGxHZDjM2AzFjU2PGZsWMTY07NiKyHWZsBmLGpsaMzYoZmxp3bERk\nO8zYDMSMTY0ZmxUzNjXu2IjIdpixGYgZmxozNitmbGrcsRGR7TBjMxAzNjVmbFbM2NQ8crrHM8Xu\nv3P/m2XviMZa+LmvqK7TZwdEdQdyZEeRP/j7dqK6PySfEtUteOk+UV3EN8+K6mZPc/+oZwCoCRsq\nquu3sYeo7q9FzUV1b80PEdU98hfZ4jzkrl+K6kZvuFtUt/Pvm0R152Q/FgRGCE73+MB1n9Yd2/6j\nGxr+oiZodzGzJJX1J8/pnoJxTuYU6Z6CkZixEZHtaF3YIoMH6BzeWD0imCWpDGwru2y0s7bdO+me\ngpG4YyMi22HGZiBmbGrM2KyYsalxx0ZEtsOMzUDM2NSYsVkxY1Pjjo2IbIcZm4GYsakxY7NixqbG\nHRsR2Q4zNgMxY1NjxmbFjE2NOzYish1mbAZixqbGjM2KGZuaR073ONI2UPn4iarWLvuy5g4TjVXc\n+6Co7sNJOaK62MIIUd26z59x2VdcVYsTfl7KvufukB1C+atdslNInuy3WlQ3aeHtorpnD5a77Cup\nuoh9lerfxWMWHRONF/T9MlFdesU6Ud3CmTtFdVO2r1I+nut9ArHero/UWDrpetF4Oa1mierSHl0k\nq/vnaPeLXD9V9O7YIq4fonN4Y0W4WNSaul5+TE7qiu0erXsKRuIzhYhsR+vCVnxojc7hjVVcVat7\nCkbaVnVR9xSMk5uTp3sKRuKOjYhshxmbgZixqTFjs2LGpsZnChHZDjM2AzFjU2PGZsWMTY07NiKy\nHWZsBmLGpsaMzYoZmxqfKURkO8zYDMSMTY0ZmxUzNjXu2IjIdpixGYgZmxozNitmbGoeOd1j7aP3\nuV3zxvoy0ViPLjgqqov1qRbVpQbJ6l4Pkp1CMmzQ3aK6P/7wrqiuyGeJqO6jv74qqnt2RCdRXeqk\nF0V1LXzPiOpynlgsqtte2FFUd7LDr0V1r89y//8eADw360tR3dKnZT+XtY9/7XbNHa9Octmn9Vdg\nxTbZImR3J0+t1z0FI5UcqNQ9BeOUZl/QPQUjcW9PRLajdWEL6BWsc3hjtW0zUPcUjNQhrJXuKRgn\nPMEjadL/O9yxEZHtMGMzEDM2NWZsVszY1LhjIyLbYcZmIGZsaszYrJixqXHHRkS2w4zNQMzY1Jix\nWTFjU+OOjYhshxmbgZixqTFjs2LGpsYdGxHZDjM2AzFjU2PGZsWMTc0j+9gxL7ZRPl5YXonOn6v7\nJs7KEo1VdH6aqC5wzXuiuq4nZf/Z7vLZ67KvwPsHdPFpr+zzWfqaaLxeny4X1U088LaobmTWaFHd\nS0+7PuVh3eYdGNQvTtn31r4A0Xj5gbLTWQYvqxDVfTnhlGy8f4tVPn6iugJHW7r+uz+/tqtovJrk\nR0R1Q77tLqpr/94D7hcddt2ldcfWOeg6ncMbq0tAL91TMJKrRa0pC6xnUWvKmLERke3Uu7AtXrwY\naWlpePrppy8/VlVVhblz52Ly5MmYN28eqqtlW3kAKCyvZy/ZhBVUbNM9BSOt27xD9xSMc6Jadkls\nd/UubElJSXjhhRd+9Fh6ejqCg4OxcOFCtGvXDunp6dd0gkRE7qp3YYuJiYG/v/+PHsvPz0dycjJ8\nfX0xbNgw5OfniwdnxqbGjE2NGZsVMzY1tzO2goIChIWFAQDCw8N/0sJGRHQtuH27R22t+595WVh+\n+PLuzJmrdQ667kcZW91+p7PZZwEALRJaNKq9+rTjcxaHto52q+1Uk30aAOCT0LpR7f07HLd7RMa1\ncqvd7dJ4zjzNuUsrqNiG0jP5GHLDOGX/ztObAQA3t+7nVtvp4mbHLyLvflGNah/YthMAENbrZrfa\nTquPO/4c2q5xbWeO5tydXdm+MmOr24+gwQCAwl3ZAIDOPRMa1T60bTsA4Ppe8W61gUgAwHcHdgEA\nbgrr2ag2sjc6/kzo7177Emem5typlVQcROvm/pfbdfuPbiwAAAT37+JWO/DSeMeqjgEA2vu1b1T7\nuyPFjr9vSIRbbSAEAFD8bSkAIKJPeKParnjVNrBSlZWVYc6cOZg/fz4AYP78+bj77rvRuXNn7N+/\nH5999tmPXlyoKyMjA+tnqe8Ru3LBq+udWavqnbgrRWnC+9iWyO5j+/Vp2X1szX/t+tOfCiq2ubwc\n7VEZJhrv9k//R1Q38UC3hr9I4a1XZPexnVpk7/vYHp3QWzbeA/cqHz9RXVHv5WjgWtmnml1IzhTV\n3RMnu48t5yX3/x3GHL4Tw4cPV/a5fSkaFRWFzMxMnDt3DpmZmYiOln+uITM2NWZsaszYrJixqdW7\nsC1YsADTp0/HwYMHMXnyZGRmZmLEiBE4evQopk6diuPHj2PEiBGemisRUaPUm7FNnTpV+fi0abLL\nvbrquxRtyuq7FG3K6rsUbaoauhRtqvjOAyKyHb5X1EDcralxt2bF3ZqaR0732PPzArdrfvbq70Rj\nldwWI6qbXPiGqO5Q0vOiunEp/xLVtV+gPvWjIUc+kp3uceibl0R1r30he5U5NUj9KldD3rjzjKiu\n/MUyUV1Q3klR3fMPthTVTZmwTFQ3K2ueqC60zd9EdZ9tullUt6J4its1Y5q77tO6YyvL4fvcVLYV\nFemegpG21K7WPQXjbCza3PAXNUHM2IjIdrQubKHdmQ+o9OrUSfcUjNTXa6juKRinf6d+uqdgJO7Y\niMh2mLEZiBmbGjM2K2ZsatyxEZHtMGMzEDM2NWZsVszY1LhjIyLbYcZmIGZsaszYrJixqXHHRkS2\nw4zNQMzY1JixWTFjU+OOjYhshxmbgZixqTFjs2LGpuaR0z2eylmgfDyrMAu9oT4D/uC0d0VjbXhc\ndrb/isoOorqBd30lqusX6foyvPrsIfSLVPcV/ZfsP/e+HbGiurlPLhTVvTBedrLyh3sPuOz7fvNR\nFPRT93dY6SMa76GwX4jqYp84Jar7j2++FNXtT1F/dkFZ9hl8kuCv7AMAPDRXNN6oebLPrIifO0pU\n1/79F90vmvSpyy6tO7benWUfbGF3Q2Ou1z0FI93Yr7/uKRgntL5FrQljxkZEtqN1YcsqzNI5vLFW\n7zmkewpG+n7zxoa/qIkpy5YdsGl33LERke0wYzMQMzY1ZmxWzNjUuGMjItthxmYgZmxqzNismLGp\nccdGRLbDjM1AzNjUmLFZMWNT446NiGyHGZuBmLGpMWOzYsamxh0bEdkOMzYDMWNTY8ZmxYxNzSOn\ne+zpvcvtmmZpmaKx9v+jpahu9Gt3iur+Z8I6UV3Bs+GiurUPDBPVffFP2eke9926VVTXc2eYqO7+\nL1yf2FCfu/bK/t0n7k4V1bWfNFlU91TkVFHduj5ForqeR2pEddVflIvqgiMeEtV17LFEVOeK1h3b\nnl3f6RzeWFuySnVPwUiVe/J1T8E4Z/fKFi67Y8ZGRLajdWGL6XmTzuGN1be37DLV7lrFROmegnFa\ndJMdsGl33LERke0wYzMQMzY1ZmxWzNjUuGMjItthxmYgZmxqzNismLGpccdGRLbDjM1AzNjUmLFZ\nMWNT446NiGyHGZuBmLGpMWOzYsamxh0bEdkOMzYDMWNTY8ZmxYxNzSOne/h88K768cMn4bNzvbJv\nc9+7RGPNmz5WVPfmrLdEdSHNN4vqwj+IcNkXnF+J8NOnlH2PHL5ONN4T/5wpqhuX94Wo7rn0g6K6\n+0tedtlXdKQSndq0Uvd99aZovD/7dhPVtbl/lKjO13uNqC7u+xPKx48dXoP2fkNc1u344xzReB2e\nPyCqe+/Dx0R1D70Y73bN2b6u+/RmbNe11Tm8sQZH8Tw2lU6R6kWtKWvfyfWi1pQxYyMi29GbsR0+\nqXN4Y63N52ceqBTtr9Q9BeMcK5Jd2todd2xEZDvM2AzEjE2NGZsVMzY17tiIyHaYsRmIGZsaMzYr\nZmxq3LERke0wYzMQMzY1ZmxWzNjUuGMjItthxmYgZmxqzNismLGpccdGRLbDjM1AzNjUmLFZMWNT\n88jpHg+3dv9Ayff/9ploLH//OFHdpKmJorrMmfeL6o5NnyCqm3VHgaju2eqJorodpT8T1e1+b4yo\n7tT9F0V1hV/LTi8Z84/Rorp3bv9UVJfx80WiutldZHuQxNrdorrauCxR3fgq2Sk5T63Kc7tm9lOu\n+7Tu2GoPyI62sbv1Ned1T8FIFTlVuqdgnLJ8ZmwqzNiIyHbqvRRdvHgxsrOzERAQgPnz5wMAPv74\nY2RkZCAgIAAAMH78eMTHu39IHAB4hd0gqrO7gT6+uqdgpIDufrqnYJzQKGZsKvUubElJSRg1ahQW\nLfpxLpCamorU1NRrOjEiIql6L0VjYmLg7+9veby2tvaqDM6MTY0ZmxozNitmbGqiV0VXrFiBTZs2\noU+fPhg5ciT8/HiJQETmcHthGzFiBO69915UVVVh+fLlWLlyJcaOrf8DVGoPHLycpzl3aV5hN8Ar\n7IYfta/sd8o977jbPNa3VaPatdW5ju/XMtatNtAOAHDgiOOTkMJCohrVXvXtLgBAUp+ebrV7XhrV\nuTtz5mp1d2t1+384vA4A0PG6QW61I5EAAPi+YCMA4MYu/RvVLt7k+BCRiFsC3Wo7HT/iuD2lXUiX\nRrWduzJnnnZlO6C7n8v+S389nNvkeDdL81vaNqp9MG8/AOCG6Ei32k7nN60FAPjeMrhR7a3HHfWJ\n7SLdajs5d2hXZmtl+Wsut+v2H805AgAI7h7iVtsp+4jj+ZAQ0r9R7e3bfgAAxPfq6FbbqaLGMX6A\nT0ij2q541TZwXVlWVoY5c+ZcfvHgSkVFRfjTn/6EmTNd30OUkZGBUR99UO8kVN7/m/rTqxryM/8Z\norpJibK3MS2eKfsQ32N9ZfexvSa8j61ZvOztSGfuE97H9oz0PrY/iOoKE2Q3NUvvY/vi9l+I6r6R\n3sd2Wz03bdXjbMS/i+qmffS6qO7cVNl9bNNGCe5j+/vPMXz4cGWf27d7lJeXAwBqamqwfv169OrV\ny+0JOTFjU2PGpsaMzYoZm1q9l6ILFixAbm4uKioqMHnyZIwbNw65ubkoKipCs2bNEBMTgxEjRnhq\nrkREjVLvwjZ16lTLY8nJyVdtcN7Hpsb72NR4H5sV72NT4zsPiMh2+F5RAzFjU2PGZsWMTc0jp3v4\njklRPl6zOwc+Pbor+9Y/5P6JIADQfMI8UV35N5NFda+kfiOqu29yhcu+U8W1OB7hpex7LLFYNF5w\n2s9FdR1/+LOorln0raK6L7ofcdmXfaYGCd19lH1fPZwtGu+ZCc1FdWeaTxPVPfFhkqhu9/F05eMX\nNu1Cs1uqXdbNLTgmG+836v+XDXm1xRuiusf6yE4hcUXrjs3VotbU9XWxqDV1CX3Vi1pT1uyWng1/\nURPEjI2IbEfrwlazO0fn8MbaUnx13otrN9lbanRPwTgXNu3SPQUjccdGRLbDjM1AzNjUmLFZMWNT\n446NiGyHGZuBmLGpMWOzYsamxh0bEdkOMzYDMWNTY8ZmxYxNjTs2IrIdZmwGYsamxozNihmbGnds\nRGQ7zNgMxIxNjRmbFTM2NY+c7pGz8jG3a5pvkh3H8vT4QFHdX7okiepWJ8s+S2DJH2Vn33f69H5R\n3Tvzt4rqRl1cLqq77bn3RHU9AmVn+59bN1FUd/aB2aK64cdk8/xq4B5R3cNehaK6qe3niOrOhQ0S\n1Z0f+TdR3fHkGLdrttfTp3XHtqn0gs7hjVV9cp/uKRhp7f5TuqdgnJJqPldUmLERke1oXdhuCffI\nlfD/Oy3bdtU9BSMNjmyjewrG6dCSzxUV7tiIyHaYsRmIGZsaMzYrZmxq3LERke0wYzMQMzY1ZmxW\nzNjUuGMjItthxmYgZmxqzNismLGpccdGRLbDjM1AzNjUmLFZMWNT446NiGyHGZuBmLGpMWOzYsam\n5pFrwa/2DlY+vq/8KMpPBSv7NvuvFY11z9uTRHUDVm4S1b32h+9FdYv7J7vsyykuR/eIDso+76HP\niMa7I2inqG7V82+J6qYnLxPV1fZ3farEdy1LcKiX+ucyvUOFaLzrB/5LVPfbT34rqiu7Tz3/hnx5\n6yzl4xd3nIZ3XKTLurEd3xeN9+r210V1fYruFNWFPej+kVRT97s+qFbrjq1rkHpRa+q6R/CMLZWb\nXCxqTZl33M26p2AkZmxEZDtaF7Z95Ud1Dm+snGKeY6/y3bYS3VMwzsUdsojB7rhjIyLbYcZmIGZs\naszYrJixqXHHRkS2w4zNQMzY1JixWTFjU+OOjYhshxmbgZixqTFjs2LGpsYdGxHZDjM2AzFjU2PG\nZsWMTY07NiKyHWZsBmLGpsaMzYoZm5pHTvfInPJ792sSxojGemHLGlHd2JyXRHUZXV8W1a1f+IKo\nLnOu7OiePa++JqoL7XJaVPeav2xxTk2rEtV999TjoroW7/xGVNf/nZWiujMfzRHVVTz2K1Fd2bM3\niuo6/mWFqK73f48W1aW9KXie3ea6S+uO7cjuDTqHN9a+U8ySVNZVHtQ9BeOsqZb94rE7ZmxEZDta\nF7aQHgN0Dm+srm2YJakManWD7ikYZ0jL1rqnYCTu2IjIdpixGYgZmxozNitmbGrcsRGR7TBjMxAz\nNjVmbFbM2NS4YyMi22HGZiBmbGrM2KyYsalxx0ZEtsOMzUDM2NSYsVkxY1Pjjo2IbIcZm4GYsakx\nY7NixqbmVVtbW3stB8jIyLiW356ImrDhw4crH7/mCxsRkacxYyMi2+HCRkS2w4WNiGyHCxsR2Y5H\nPvOgrtzcXCxZsgQ1NTVISUlBSkqKjmkYZ8qUKfDz84O3tzd8fHwwe/Zs3VPSYvHixcjOzkZAQADm\nz58PAKiqqsKbb76JwsJCREZG4oknnkDLli01z9RzVD+Tjz/+GBkZGQgICAAAjB8/HvHx8TqnaQwt\nC9vSpUumOPuUAAACUElEQVSRlpaGkJAQvPLKKxg4cODlf5ymbsaMGWjdumnfTZ6UlIRRo0Zh0aJF\nlx9LT09HcHAwnnzySSxbtgzp6ekYO3asxll6lupnAgCpqalITU3VNCtzefxStLKyEgAQGxuLkJAQ\nxMXFIT8/39PTMBbvvgFiYmLg7+//o8fy8/ORnJwMX19fDBs2rMk9Z1Q/E4DPF1c8vmPLz89HeHj4\n5XZ4eDjy8vKQkJDg6akYx8vLCy+//DJCQ0MxbNgwJCYm6p6SMQoKChAWFgbA8ZxpagubKytWrMCm\nTZvQp08fjBw5En5+frqnZAQtl6KkNnPmTAQFBaGkpARz5sxBVFQUAgMDdU/LCNyZWI0YMQL33nsv\nqqqqsHz5cqxcubJJXZ7Xx+OXolFRUSgtLb3cLikpQXR0tKenYaSgoCAAQIcOHZCYmIitW7dqnpE5\nrnzelJaWokuXLppnpF/btm3h5eWFVq1aYeTIkfj22291T8kYHl/YWrVqBcDxymhZWRl27tyJqKgo\nT0/DOGfPnkVVleNT0CsqKrBjxw6+wnWFqKgoZGZm4ty5c8jMzOQvQwDl5eUAgJqaGqxfvx69evXS\nPCNzaHmvqPN2jwsXLiAlJQWjR4/29BSMU1ZWhnnz5gEA2rRpg0GDBiE5OVnzrPRYsGABcnNzcerU\nKbRt2xb33XcfbrnlliZ9u4fzZ1JRUYHAwECMGzcOubm5KCoqQrNmzRATE4N77rmnyb+i7sQ3wROR\n7fCdB0RkO1zYiMh2uLARke1wYSMi2+HCRkS2w4WNiGznfwGsKg8AmsmeSgAAAABJRU5ErkJggg==\n"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD/CAYAAADWreLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAA8xJREFUeJzt1jGrznEYx2FHDyspmZhlUkYlyxksdFZvQYwHu8Fq8Crs\nstqkFK9AYTPwDv72Zzs+/brPqeva7/pOn+6Dbdu2cwD/6fz0AOBsExEgEREgEREgEREgEREgEREg\nEREgEREgEREg2U0P2Pfkw8H0hKV+Prs2PWGpH3ceTk9Y6vGVB9MTljp+e3TiG58IkIgIkIgIkIgI\nkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgI\nkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgI\nkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkOymB+z7fPPi9ISljg+fT09Y6uWn\ne9MTlnp6/df0hFPHJwIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIku+kB+95duDQ9YalbH4+mJyx1+f7V6QlLvTl8NT1hqRfnHp34xicCJCICJCICJCICJCIC\nJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCIC\nJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCIC\nJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJLvpAfu+/v0zPWGp99++T09Y6vWX\nu9MTlvp94/b0hFPHJwIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkB9u2bdMjgLPLJwIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIk/wBVaCFzoNoUHAAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 203
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Train the map"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import time\n",
      "start = time.time()\n",
      "for i in range(num_iterations):\n",
      "    \n",
      "    # Pre-calculate the number of iterations (which will never be so impossibly large as to not store the indices)\n",
      "    idx = np.random.randint(0, n_samples, (num_iterations,))\n",
      "    \n",
      "    # Pick a random vector \n",
      "    this_sample = X[idx[i],:]\n",
      "\n",
      "    # Figure out who's the closest weight vector (and calculate distances between weights and the sample)\n",
      "    find_winner(som_weights, this_sample, abs_weight_distances)\n",
      "    winner_idx = np.argmin(abs_weight_distances)\n",
      "\n",
      "    # Calculate the new learning rate and new learning spread\n",
      "    normalized_progress = float(i) / float(num_iterations)\n",
      "    learning_rate = scale(learning_rate_initial, learning_rate_final, normalized_progress)\n",
      "    learning_spread = scale(learning_spread_initial, learning_spread_final, normalized_progress)\n",
      "    \n",
      "    # Update those weights\n",
      "    update_weights(this_sample, som_weights, winner_idx, grid_indices, learning_rate, learning_spread)\n",
      "    \n",
      "print \"Training %d iterations on a %d-square grid took %f seconds\" % (num_iterations, grid_size[0], time.time() - start)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Training 200 iterations on a 20-square grid took 0.033242 seconds\n"
       ]
      }
     ],
     "prompt_number": 204
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Show off the finished Kohonen map"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "the_som = np.reshape(som_weights, (grid_size[0], grid_size[1], n_dims))\n",
      "figure(figsize=(5,5))\n",
      "imshow(the_som);\n",
      "\n",
      "figure(figsize=(4,4))\n",
      "axis('off')\n",
      "numRows = 4\n",
      "imshow(X[:,np.newaxis,:].reshape(numRows,X.shape[0]/numRows,3))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 205,
       "text": [
        "<matplotlib.image.AxesImage at 0x115ead390>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAEzCAYAAAC7cS8aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHw1JREFUeJzt3VtsVNe5B/Bv75nx2MZXgg32cEmCaWPShktMLk3PERcJ\n4sjqQ9P06OSpLzygqEJVj3qkVJFQaBURhCpadHhAraJWRzpSHvqaEllIPeVwMyQGbBMYxw54bHwB\n4/F47jP7PBgQ9v722PPFzFpd/v/e1mwv1vLMsGbt/2x/23IcxyEAAIPYqicAALDUsLABgHGwsAGA\ncbCwAYBxsLABgHGwsAGAcfxPe4DOzs6nPQQALFN79uxhH3/qCxsR0X98dL/oPol4g2isdLZc1C+f\n94n6ZfIBUb90Pijql3VkL5llW6J+FcG8qN/K6rSo34bauKjfs9XTon6rgklRv0orK+pXTrLnsy4r\nu9y0KSMbr0n28lFdVvY+K8sX3+/agVrPY0pPRWOTPSqH11YmdlH1FLR0b/yy6ilop2/kuuopaEm8\nY+vt7aVTp05RLpej9vZ2am9vX8p5AQCIiXdsn3zyCe3fv58++OAD+tvf/kbRaLTof6Oq/kXp8EYL\nVL2iegpaeqbhZdVT0E5r0/dUT0FLooUtHp/NQTZv3kwNDQ20ZcsWCofDSzoxAAAp0cIWDocpFAo9\nbodCIbp161bR/w4yNh4yNh4yNjdkbLySfHnw5AIWm+wpqh2f+pLiU18uup2IXqFE9Iq4nZy+TMnp\ny4tup6a7KDXdJW5nYhfnLGSZ2EXKJm54Hs/GLlA2dkHczsQuUKaIdip6gVJReXtm8hLNTF5adPve\n+OU5C1ix7cjdborc7V50e3D4Gg0OXxO3+4evU//w9UW3bw1fp1vfot03cn3OYnb73sCc9vzj3Xd7\nqPtuj7h9abSHLo0uvn1urIfOjcnbZ8ev09nx64tue7EkZYvi8TgdOnSIPv74YyIi+tOf/kRbt26l\n7du3u362s7MTl3swcLkHD5d78HC5h9u1A7We17GJdmyVlZVENPvN6NjYGF29epVaWlok/xQAwJIT\nX+7xs5/9jE6dOkXZbJba29uppqam6H8jNtmDb0YZmdhFfDPKuDd+Gd+MztM3ch3fjDLEC9vmzZvp\nd7/73VLOBQBgSSj9ywPs1njYrfGwW3PDbo1Xkr8VDTYUHwjHJ2RhsJ2UpZ5OVhbm2znZlxV2TvZl\nhe3IPossn2w8KpON5wRk42X9srdkxpKNlyNZ2C2L5ImkNxhxLNmIGTsn6pf0yWYq++qHKJ+TvQ5e\nlO7YomM3Fv6hZSg9fWnhH1qG7o91LfxDy8yNYVwLykE9NgAwjtKFrabxBZXDa6useofqKWhpZWOb\n6ilo54Vm5NQc7NgAwDjI2DSEjI2HjM0NGRsPOzYAMA4yNg0hY+MhY3NDxsbDjg0AjIOMTUPI2HjI\n2NyQsfGwYwMA4yBj0xAyNh4yNjdkbDzs2ADAOMjYNISMjYeMzQ0ZG68k1T3syhH2cSt4z/NYeU1M\nNFY6UCHqZ6VXiPo56TpRP1/Gu5qBHUiRryzBjyf8LLLKhNVEKoVVMypl84wHvKtKJP1lFA/w1VSm\n7YxovGoSlvi2ZFUzLEc2Hlkev5+V9T5GRBmfbLwYyX6/vHCrFPQv7R5L6Y6ttul5lcNrK1iDumOc\nutWoUzffC82bVU9BS8jYAMA4She2qZGvVQ6vrVQU98/kPBjF/VbnuzHcq3oKWsKODQCMg4xNQ8jY\neMjY3JCx8bBjAwDjIGPTEDI2HjI2N2RsPOzYAMA4yNg0hIyNh4zNDRkbDzs2ADAOMjYNIWPjIWNz\nQ8bGw44NAIyDjE1DyNh4yNjckLHxSlLdw994u/g+dqVssBm++sNC7FStbLxUUtTNScmqZpBVJutX\nLqt6Qitk46WDss/MqF/2vDyw86J+NZQW9atyZK87CatmkC2bZ9aS9ZuxZfNM+7yrsxQScAyq7vHg\ndkTl8NpKTnarnoKWYiPnVU9BOzeGUdOQg4wNAIyjdGGrWx9SOby2yuu3qJ6ClqqaXlM9Be280Iz7\nhnCwYwMA4yBj0xAyNh4yNjdkbDzs2ADAOMjYNISMjYeMzQ0ZGw87NgAwDjI2DSFj4yFjc0PGxsOO\nDQCMg4xNQ8jYeMjY3JCx8bBjAwDjIGPTEDI2HjI2N2RsvJJU96BnhvjHH0wQPcNXA/Bna2Rj+WRV\nQexMXNYvLasqYaW9q23kslEKrrzP93OqRePlgpaoHwUDom5ZYZWOuOX9WZv0l5Ed4Ku3TDkZ0XjT\nluz3W1lgnoU4VlbUz3K8qnRkyCpQwSMvrO6RJtn7OmPLqntI3i2F/qerzdg2rlI5vLYqG19UPQUt\nVTb/QPUUtNMa+q7qKWgJGRsAGEdtxtY/oXJ4bcXHelRPQUvx4f9TPQXt9EW+Uj0FLWHHBgDGQcam\nIWRsPGRsbsjYeNixAYBxkLFpCBkbDxmbGzI2HnZsAGAcZGwaQsbGQ8bmhoyNhx0bABgHGZuGkLHx\nkLG5IWPjYccGAMZBxqYhZGw8ZGxuyNh4Janu4VszXHSfgG9GNtakrPpFPi0bL5DNifrZGe/qHoVY\n6TpRv5RTL+qXE370OZasmkhO+FmbEtWHIEo4sv8CGZJVsXBsWXUPx5G9z6T9csLqHsKXnfKO7Pn0\nojZj65tWOby2YkO3VE9BS4nIWdVT0M6NyE3VU9ASMjYAMI74VPS9996jiooKsm2bfD4fffTRR0X/\nG3WtstNG01Wt3aR6ClqqCL2hegraeSH0HdVT0NK3ytgOHTpEVVVVSzUXAIAl8a1ORZ1vGfghY+Mh\nY+MhY3NDxsYT79gsy6IPP/yQGhsbadeuXdTW1raU8wIAEBMvbIcPH6b6+noaGhqiI0eOUEtLC9XV\nFXcpAjI2HjI2HjI2N2RsPPGpaH397HVRa9eupba2Nurq6vL82SdPOR/0TRfVngpP0lR4cvHtwQma\nGpwQt6O3Ryl6e3Tx7TtDFL0zJG7Hhm7NOfVcqB0f7qb4cLe4nRy+QMnhC4tupyJnKfXEKWCx7WTk\nH5SM/GPR7UTk7JxTzmLb05HzNB05v+j2eOQyjUcui9t3ItfoTuTaotvhSB+FI33i9o3IzTmnnwu1\neyK3qCdy61u0w9QTCZesfT0SputFtL1YjiAoS6VSlM/nqaKigqLRKB06dIjef/99WrXK/ZcEnZ2d\n9JuR/2T/nQd90567tsx4bbHTIiKivPgC3RWyflnZX0+kMhs9j8WGbnnu2pJp2Z2/U86zon45W3ZB\nsGXJTgYC5H1BaSJy1nPXVk8x0XjrnXuifhtpdOEfYoR8svHqHf72kDciNwvu2io8b9tXmC29QFfU\ni8gS5PW1b79Ee/bsYY+J3n1TU1N09OhRIiKqrq6mjo4OdlEDAFBBtLA1NjY+Xti+DWRsPGRsPGRs\nbsjYePjLAwAwDv5WVEO4jo2H69jccB0bryTVPezmcfZxayxJdnOSPRaskoXBVC/7EsBKVIr65dOy\nxbks5f2Zkk/epeoQ/9LYMzWi8ay47MuYdCYo6uc4ZaJ+Psc7tLbtLPl8fBjuWCnReJadEfXLZ/n3\n7UIyPlm1Dcp7zNPKEhX4HZy8tCqIqJuw5gmRZRtU3aN+a7nK4bVVu3Gt6iloqXLt66qnoJ0XkMey\nkLEBgHGULmyTX8q286ab6h9a+IeWofjQOdVT0M4N5LEs7NgAwDjI2DSEjI2HjM0NGRsPOzYAMA4y\nNg0hY+MhY3NDxsbDjg0AjIOMTUPI2HjI2NyQsfGwYwMA4yBj0xAyNh4yNjdkbDzs2ADAOMjYNISM\njYeMzQ0ZG68k1T2smpmi+9jBhGisskpZtY1AXFaNgjKy8dJJ2WeKf0ZWvSQwLavSkZyRlZZ2srIq\nJIGcrBpFhVP8e4yIyBKWFI9XTon6ZdOTC/8QI+HInhe/7OWjgCMr8m0Ji4Nnhb+fF7UZ22VZyRjT\n3eubWPiHlqHp2943DFqubt7pVz0FLSFjAwDjqM3YXg6oHF5bz7Tixjic6vW4Kfd831nnfbez5Qw7\nNgAwDjI2DSFj4yFjc0PGxsOODQCMg4xNQ8jYeMjY3JCx8bBjAwDjIGPTEDI2HjI2N2RsPOzYAMA4\nyNg0hIyNh4zNDRkbDzs2ADAOMjYNIWPjIWNzQ8bGK011D48SA5YvS5bfYY8F+IcXVJ6TFa8M2j5R\nPzubFfXLrPCutpGoiVL1M/yiH6ytFI1XmZQ9odmEbJHNJutE/fwp72N2bJBWNlWwx4KWrIxFeVbY\nLzMm6hfPy16H8jj/vk4FpylRed+7o/B9XZ4TVp/JysbLScuQeFCbsbWVZF39p9P4kqzkj+lWbmpV\nPQXttDy7XvUUtISMDQCMozZj65Kdxplu7GpU9RS0dP9Wn+opaCc8eFv1FLSEHRsAGAcZm4aQsfGQ\nsbkhY+NhxwYAxkHGpiFkbDxkbG7I2HjYsQGAcZCxaQgZGw8ZmxsyNh52bABgHGRsGkLGxkPG5oaM\njYcdGwAYBxmbhpCx8ZCxuSFj45VkZVmRK/6Usywnq4JQ4eRF/QLCKghWUNavzB8T9avwj4v6VZMl\n6pfP3hP1o5SsCkkuVSbqF4jJ+gXjstJZudyMqN90QlqqS/Y6VKdlz0ss6119phDbJ/v/l8kliu6z\nlvZ4z0M0iyUycTmncnhtjXwh+09jutFryJPm++ZrWfkk0yFjAwDjKF3YVr0sO40zXdO2FaqnoKXV\n30eeNN+G5xtVT0FL2LEBgHGQsWkIGRsPGZsbMjYedmwAYBxkbBpCxsZDxuaGjI2HHRsAGAcZm4aQ\nsfGQsbkhY+NhxwYAxkHGpiFkbDxkbG7I2HjYsQGAcZCxaQgZGw8ZmxsyNl5Jqns0ZPkFLJvLU4NH\n4Q87K6sS4JMVBSHbklW/sAKyAfPl3r9fWdCh8gqP4+Vp0XhOYErUz7JSon6UEcYMmXLPQ/GRe1Qf\nCrDH8jOyt7ITE75+40lRv0RW1i8T56t7RJMP6H7c+3eICt/XvoRs0+ET1o7NWktbdFbpjq1pO86E\nOaGXZSV/TBdqW6N6CtoJPVunegpawsoCAMYpuLCdPHmS9u/fT7/85S8fP5ZIJOjjjz+mAwcO0NGj\nRymZlG2tiYhGrshON00XuRxXPQUtRbruqp6CdiKDD1RPQUsFF7adO3fS+++/P+ex06dP06pVq+j3\nv/89rVy5kk6fPv1UJwgAUKyCC1trayutWDH3mqpwOEy7d++mQCBAu3btonA4LB4cGRsPGRsPGZsb\nMjZe0StLf38/NTc3ExFRKBT6VgsbAMDTUPR35I5T/NfjI1fyj3dnj3K1pu32nIxt/vHQ92cfH/5i\ndrzmbdai2pEvZ9uhrcW1122hh+2H42+1F9f+YvbmHKFtgaLaTa8/HO9hnvZolxa5HKeJm0na8u8r\n+eOXZjOV0I66otrNP6iebV+cvWwg9Mozi2tfmL15TOjVhuLa29c8HH/s4XwaF9d+mKM92p092X4y\nY5t/vKl1LRERDX8xMvv7bmtaXPvqw/ZLxbXXN9XPjt83Ojuf1tWLao/cnH1+mr7TUFR7Q9XD99vD\nTO3RTq37/BCtWlP1uD3/+NDt2fba9cW1N9TMvl/uDM9eJrSuuXZx7ZGH7abi2k3Ns2eGQ8Oz99Rd\n21yzqLYXy1lgpRobG6MjR47QsWPHiIjo2LFj9OMf/5iee+45+vrrr+mvf/3rnC8X5uvs7KT/qXqL\nPfbkgjefnRFexya8HMZ2hHepCsj+/Clf7v2nMJHLce/T0QL9CnEC1aJ+luV9XVlBT+E6tkjXXc/T\n0VJfx+YXXseWisn6lQ3z17FFBh8UPB31Z4TXsT3Q/zq2jn/7De3Zw9+pquhT0ZaWFjpz5gyl02k6\nc+YMbdq0qegJPYKMjYeMjYeMzQ0ZG6/gynL8+HH64IMPaGRkhA4cOEBnzpyhvXv30sTEBB08eJDu\n379Pe/fuLdVcAQAWpeD+/eDBg+zjv/rVr5Zk8EKnostZwVPRZazQqehytdCp6HKFVQUAjIO/FdUQ\ndms87NbcsFvjlaS6xzqr+G84M8I1Lyv8Mi5vyQbMl/HVJhbsVy576vMVsm+5yC8sESWt7hGUPZ9W\nQva1WjYve+GzWdm3oo4tq7Lii8n+XC6Xkr0OM/dion7OpKgb5ZOyeeaFVXm8KN0yDVxe4t/GEHcu\noR4b587FcdVT0M7w7ajqKWgJ54IAYBylC9tzLwtPqwy3bgfuecBZ90qD6ilop3l94Svwlyvs2ADA\nOMjYNISMjYeMzQ0ZGw87NgAwDjI2DSFj4yFjc0PGxsOODQCMg4xNQ8jYeMjY3JCx8bBjAwDjIGPT\nEDI2HjI2N2RsPOzYAMA4yNg0hIyNh4zNDRkbryTVPUJl/OOxgPextLBKRyovW6tTwqciWebxCyw0\nXoF+AX+Kgh7HM37ZE+P4ZKf9Dsk+fKy88GbYAe/qHrY/R7bH8TKfrHpJQFBrn4goZ8uqWCSzsqog\n6Rj/YZdMJGkm5v2eSI/JxstPCat0JGSvuyN9v3hQumP7bhsyNs6GV2U3XjHduldXqZ6CdlaHkMdy\nkLEBgHGULmxfdSFj43xzYVr1FLR058KE6iloZzSCPJaDHRsAGAcZm4aQsfGQsbkhY+NhxwYAxkHG\npiFkbDxkbG7I2HjYsQGAcZCxaQgZGw8ZmxsyNh52bABgHGRsGkLGxkPG5oaMjYcdGwAYBxmbhpCx\n8ZCxuSFj45WkukeNv/gFLCes7pFxZGt10pI9FTN+WXWPmE/WLyGcZ1b4GeY4wg8fYbEGKyuLJ+yM\ncLy07PfLpAKifr5UQtQvH5M9L7lpWdUTab98WvjC55c2llK6Y+u5tLSlSkwxcAE1tji3LyJjm2/0\nLjI2DjI2ADCO0oXtxR1YVznPvYo69pz1ryBjm2/1GmRsHKwsAGAcZGwaQsbGQ8bmhoyNhx0bABgH\nGZuGkLHxkLG5IWPjYWUBAOMgY9MQMjYeMjY3ZGw87NgAwDjI2DSEjI2HjM0NGRsPKwsAGAcZm4aQ\nsfGQsbkhY+OVpLpHwOYrKPgty/uYsKpEwJKVBQnYQVE/n192KmD5vPtV2mmq9jhuU4VovGRe9hmW\nywmrLmRk/ayEdz9/yqJAgn9f2B6PLygte17yJCsnIih0Q0REdpb/r2rnfJ7HiIgcYbWUvPD1k45H\njkHVPb73Cs6EOS2v1auegpY2tCFjm2/16irVU9ASVhYAMI7She36RWRsnPD5SdVT0NI3XcjY5hsd\njamegpawYwMA4yBj0xAyNh4yNjdkbDysLABgHGRsGkLGxkPG5oaMjYcdGwAYBxmbhpCx8ZCxuSFj\n42FlAQDjIGPTEDI2HjI2N2RsPOzYAMA4yNg0hIyNh4zNDRkbryTVPcrs4hcwR7jmOpasSkfALyvu\n6A80yPr5Vov6BRxZdY+ZbE7UL5VOifrlE7LqFxSXzdNKlbYaRd4n+68T9AVE/eI+2f8Hi4TlRJa2\n2EbJKd0yXUXGxuo7N6p6CloauDKuegrauXt3WvUUtIRzQQAwTsH99MmTJ+nKlStUU1NDx44dIyKi\nTz/9lDo7O6mmZvbU7d1336WtW7eKBn8JGRur9XXZaarpntsuO+032Zo11aqnoKWCC9vOnTvpzTff\npBMnTsx5vKOjgzo6Op7qxAAApApumVpbW2nFCneJameJyvgiY+MhY+MhY3NDxsYTfbXz2Wef0fnz\n52nHjh20b98+qqiQfVMHAPA0FB1y7d27l06cOEG//vWvaXR0lD7//PMF+zy5M7t6Mf+4/dIr9pz2\n/ONERNcu5ujaxVwR7Qxdu5iRty/E6dqF+KLbPecmqefcpLjdd250zg5t/m5t/vGb50bo5rkRcbv/\n/F3qP3930e2Bi6M0cHFU3u4ap4Gu8cW3r4zP2Zk92X5ue0PB40REA90TNNA9sfj2tQkauCZvf9N3\nj77pu7fo9tDAAxoaeCBu37077dqlPdmef3x8KkHjUwlxe2ImRRMzqcW342maiKdL1vZiOQucV46N\njdGRI0cef3nwpMHBQfrjH/9Ihw8f9uzf2dlJqVDxeZz8OrZyUb+8X3bxZzqwTtQvLryOLVbq69hS\n/yTXscWF17HFZXFIOiqb58w3stvlTX41JeoXG5D1y0wmRf2crDBeEsRbB//rv2nPnj3ssaJXj8nJ\n2d1GLpejs2fP0rZt24qe0CPI2HjI2HjI2NyQsfEKZmzHjx+n3t5eikajdODAAXrnnXeot7eXBgcH\nye/3U2trK+3du7dUcwUAWJSCC9vBgwddj+3evXvJBsd1bDxcx8bDdWxuuI6Nh5UFAIyDvxXVEDI2\nHjI2N2RsvJJU96jyl7GPV/iyVOXnp+AIp5aza2X9yppE/QJlz8r62d6nm1UBi2qDG/h+OVl1iLJc\nYuEfYsSyUVG/VEpWADGf8v42zp/xUSDlY485admHpCP8bPVZso52kJ//gv3K+H5WwOd5jIjI9gur\ngtiyqiDyS/eFVUg8KN2xbXu1JOvqP53vvcEvasvd81uRsc3X1ISMjYOMDQCMo3Rh++JCVuXw2rp+\n9hvVU9DS118iY5tvZAQZGwc7NgAwDjI2DSFj4yFjc0PGxsOODQCMg4xNQ8jYeMjY3JCx8bBjAwDj\nIGPTEDI2HjI2N2RsPOzYAMA4yNg0hIyNh4zNDRkbDzs2ADAOMjYNIWPjIWNzQ8bGK8nKUuNR3aOQ\nnPDeBRnfSlm/4HpZv7LnRf18tuweC76crAqCnZOdsji2sDpEXnZPgHROFk/k8rJ7LOQd2fNpCfcE\nliWr7uELyKq6ULD4/3tERBSQPZ9WTnjvifzS3NLzEaU7tq7zwht+GK77f79SPQUt9XePqZ6CdkaG\nZTdrMR0yNgAwjtKFre014fbacFv+5buqp6CljVsaVU9BO03NssKqpsOODQCMg4xNQ8jYeMjY3JCx\n8bBjAwDjIGPTEDI2HjI2N2RsPOzYAMA4yNg0hIyNh4zNDRkbDzs2ADAOMjYNIWPjIWNzQ8bGw44N\nAIyDjE1DyNh4yNjckLHxSlLdI+jjTzkDdt7zmLS6hxWoF/XLB5qE/daI+pHtfQrh89eSP/CMRz9Z\nFYRcmeylzvhTsn52XNQvS0nPYzb5yCaP94sjqyaSJ2FVCWE1CseWVRPJl/NVOvJlAc9jRETWiqBo\nPDsmez4zjqx6CeXysn4elO7YXnld9qSbbuu/vqh6ClrauEX4IWKwprV1qqegJWRsAGAcpQvbxXOy\n0xzTffn3HtVT0FJ/913VU9DOyNAD1VPQEnZsAGAcZGwaQsbGQ8bmhoyNhx0bABgHGZuGkLHxkLG5\nIWPjYccGAMZBxqYhZGw8ZGxuyNh42LEBgHGQsWkIGRsPGZsbMjYedmwAYBxkbBpCxsZDxuaGjI1X\nkuoeluVdfcCzj69SNphP+EJL+9myeUqeEyIi2yerKhHwyT5EfH5ZlRW/8PXzU0zUL2PJqlE4WVlV\nCSclez4zwmoiAUdW4stXXiHqFysXPi9+2XiBpOz186I4Y/MuUbOcffn3a6qnoKVw94jqKWhneGhS\n9RS0hIwNAIyjOGOTneaYbuu/fl/1FLTUskVWDNRkzWtlhVVNhx0bABgHGZuGkLHxkLG5IWPjYccG\nAMZBxqYhZGw8ZGxuyNh42LEBgHGQsWkIGRsPGZsbMjYedmwAYBxkbBpCxsZDxuaGjI2HHRsAGAcZ\nm4aQsfGQsbkhY+NZjuPIyg0sUmdn59P85wFgGduzZw/7+FNf2AAASg0ZGwAYBwsbABgHCxsAGAcL\nGwAYpyT3PJivt7eXTp06Rblcjtrb26m9vV3FNLTz3nvvUUVFBdm2TT6fjz766CPVU1Li5MmTdOXK\nFaqpqaFjx44REVEikaA//OEPNDAwQM8//zz9/Oc/p/Ly5XOBN/ecfPrpp9TZ2Uk1NTVERPTuu+/S\n1q1bVU5TG0oWtk8++YT2799PDQ0N9Nvf/pbeeOONxy/Ocnfo0CGqqqpSPQ2ldu7cSW+++SadOHHi\n8WOnT5+mVatW0S9+8Qv685//TKdPn6Yf/ehHCmdZWtxzQkTU0dFBHR0dimalr5KfisbjcSIi2rx5\nMzU0NNCWLVsoHA6XehrawtU3RK2trbRixYo5j4XDYdq9ezcFAgHatWvXsnvPcM8JEd4vXkq+YwuH\nwxQKhR63Q6EQ3bp1i7Zv317qqWjHsiz68MMPqbGxkXbt2kVtbW2qp6SN/v5+am5uJqLZ98xyW9i8\nfPbZZ3T+/HnasWMH7du3jyoqZLe/M42SU1HgHT58mOrr62loaIiOHDlCLS0tVFeHG+ISYWfC2bt3\nL/3kJz+hRCJBf/nLX+jzzz9fVqfnhZT8VLSlpYUikcjj9tDQEG3atKnU09BSff1spYa1a9dSW1sb\ndXV1KZ6RPp5830QiEdq4caPiGalXW1tLlmVRZWUl7du3jy5duqR6Stoo+cJWWTl7h/De3l4aGxuj\nq1evUktLS6mnoZ1UKkWJRIKIiKLRKHV3d+Mbrie0tLTQmTNnKJ1O05kzZ/BhSESTk7N/AJ/L5ejs\n2bO0bds2xTPSh5K/FX10uUc2m6X29nZ66623Sj0F7YyNjdHRo0eJiKi6upp++MMf0u7duxXPSo3j\nx49Tb28vTU9PU21tLf30pz+l1157bVlf7vHoOYlGo1RXV0fvvPMO9fb20uDgIPn9fmptbaW33357\n2X+j/gj+CB4AjIO/PAAA42BhAwDjYGEDAONgYQMA42BhAwDjYGEDAOP8P48/gklDZSWIAAAAAElF\nTkSuQmCC\n"
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAREAAAD/CAYAAADWreLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAA8xJREFUeJzt1jGrznEYx2FHDyspmZhlUkYlyxksdFZvQYwHu8Fq8Crs\nstqkFK9AYTPwDv72Zzs+/brPqeva7/pOn+6Dbdu2cwD/6fz0AOBsExEgEREgEREgEREgEREgEREg\nEREgEREgEREg2U0P2Pfkw8H0hKV+Prs2PWGpH3ceTk9Y6vGVB9MTljp+e3TiG58IkIgIkIgIkIgI\nkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgI\nkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgI\nkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkOymB+z7fPPi9ISljg+fT09Y6uWn\ne9MTlnp6/df0hFPHJwIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIku+kB+95duDQ9YalbH4+mJyx1+f7V6QlLvTl8NT1hqRfnHp34xicCJCICJCICJCICJCIC\nJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCIC\nJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCIC\nJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJCICJLvpAfu+/v0zPWGp99++T09Y6vWX\nu9MTlvp94/b0hFPHJwIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkB9u2bdMjgLPLJwIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIkIgIk\nIgIkIgIkIgIkIgIkIgIkIgIk/wBVaCFzoNoUHAAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 205
    }
   ],
   "metadata": {}
  }
 ]
}