som-random.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import patches as patches\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_data = np.random.randint(0, 255, (3, 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "network_dimensions = np.array([5, 5])\n",
    "n_iterations = 10000\n",
    "init_learning_rate = 0.01\n",
    "\n",
    "normalise_data = True\n",
    "normalise_by_column = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = raw_data.shape[0]\n",
    "n = raw_data.shape[1]\n",
    "\n",
    "# initial neighbourhood radius\n",
    "init_radius = max(network_dimensions[0], network_dimensions[1]) / 2\n",
    "# radius decay parameter\n",
    "time_constant = n_iterations / np.log(init_radius)\n",
    "\n",
    "data = raw_data\n",
    "if normalise_data:\n",
    "    if normalise_by_column:\n",
    "        col_maxes = raw_data.max(axis=0)\n",
    "        data = raw_data / col_maxes[np.newaxis, :]\n",
    "    else:\n",
    "        data = raw_data / data.max()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "net = np.random.random((network_dimensions[0], network_dimensions[1], m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_bmu(t, net, m):\n",
    "    \"\"\"\n",
    "        Find the best matching unit for a given vector, t\n",
    "        Returns: bmu and bmu_idx is the index of this vector in the SOM\n",
    "    \"\"\"\n",
    "    bmu_idx = np.array([0, 0])\n",
    "    min_dist = np.iinfo(np.int).max\n",
    "    \n",
    "    # calculate the distance between each neuron and the input\n",
    "    for x in range(net.shape[0]):\n",
    "        for y in range(net.shape[1]):\n",
    "            w = net[x, y, :].reshape(m, 1)\n",
    "            sq_dist = np.sum((w - t) ** 2)\n",
    "            sq_dist = np.sqrt(sq_dist)\n",
    "            if sq_dist < min_dist:\n",
    "                min_dist = sq_dist # dist\n",
    "                bmu_idx = np.array([x, y]) # id\n",
    "    \n",
    "    bmu = net[bmu_idx[0], bmu_idx[1], :].reshape(m, 1)\n",
    "    return (bmu, bmu_idx)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def decay_radius(initial_radius, i, time_constant):\n",
    "    return initial_radius * np.exp(-i / time_constant)\n",
    "\n",
    "def decay_learning_rate(initial_learning_rate, i, n_iterations):\n",
    "    return initial_learning_rate * np.exp(-i / n_iterations)\n",
    "\n",
    "def calculate_influence(distance, radius):\n",
    "    return np.exp(-distance / (2* (radius**2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_iterations):\n",
    "    # select a training example at random\n",
    "    t = data[:, np.random.randint(0, n)].reshape(np.array([m, 1]))\n",
    "    \n",
    "    # find its Best Matching Unit\n",
    "    bmu, bmu_idx = find_bmu(t, net, m)\n",
    "    \n",
    "    # decay the SOM parameters\n",
    "    r = decay_radius(init_radius, i, time_constant)\n",
    "    l = decay_learning_rate(init_learning_rate, i, n_iterations)\n",
    "    \n",
    "    # update weight vector to move closer to input\n",
    "    # and move its neighbours in 2-D vector space closer\n",
    "    \n",
    "    for x in range(net.shape[0]):\n",
    "        for y in range(net.shape[1]):\n",
    "            w = net[x, y, :].reshape(m, 1)\n",
    "            w_dist = np.sum((np.array([x, y]) - bmu_idx) ** 2)\n",
    "            w_dist = np.sqrt(w_dist)\n",
    "            \n",
    "            if w_dist <= r:\n",
    "                # calculate the degree of influence (based on the 2-D distance)\n",
    "                influence = calculate_influence(w_dist, r)\n",
    "                \n",
    "                # new w = old w + (learning rate * influence * delta)\n",
    "                # where delta = input vector (t) - old w\n",
    "                new_w = w + (l * influence * (t - w))\n",
    "                net[x, y, :] = new_w.reshape(1, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEICAYAAABBKnGGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFR9JREFUeJzt3XmUXGWdxvHvYxJIyELUtAgkBBkQ\nDCCgTUABB8EFJYoeN1BwcJgTdVxwxtERN8TB9YwI6hw9GRbZGQbEhUEWBYSoLAmgEBIWIRDWNEuW\nZksIv/njfVuri67qzlvVfavJ8zmnT6rrLu+vbt167nvfe6ujiMDMbH29qOoCzGx0cniYWRGHh5kV\ncXiYWRGHh5kVcXiYWZFhDQ9JIWnb/HiCpF9JWinpf4ez3XaStEjSvoPMs5WkXkljRqisSknaXtJN\nklZL+nTV9VSlE953SR+SdGkljUdE0x9gb+APwErgMeD3wO6DLZeXDWDb/Pgw4DpgbJP5BXwOuAN4\nCrgX+Baw8VDaeyH+AD/N2/Gguue/n58/vIKaTgK+X1fjsW1u443AFXm/WzrA9K3z9CeBJcCb6qb/\nC/AQsAo4uXYfamXZQWq+EvinYdzuW+f3vOFnaCR/mvY8JE0BLgR+CLwE2BI4Bnim2XINzARuj4hn\nm8zzA2Au8GFgMvA2YH/g3CY1bghH+9tJ2wQASWOB9wN/qaiemcCidq0sv556T5A+uJ9rsNjZwI3A\nS4EvAedJ6srreyvwBdK+MxPYhrTftmPZYTPq9uVBkq4bWDHIPP8ILAYeBy4BZtZMC2Bb0sZfA6wF\neoEjBljPdsA6YHbd8zNIYbVfzVHux8BFpB3sTaSd4FekI8X1wLHA/Jp1nAAsy9MXAvvUTPsaKZxO\nA1aTPhTdNdOXko9MwGxgQV7Pw8BxAx0RSEeg/yD10lYDlwLTatb5YeAe4FHgK7VtDLBdfgr8Z27v\nxfm5OcCvgfnkngfwd8DleZ2PAGcCU+tex1HArfm9OgUY36DNhuvKz68Dns7v5dz8vq7Jv/8qz7cF\ncD7QA9wNfLpum58HnJG3ZcOjdX5/l9Y998q8T0yuee5q4GP58VnAN2um7Q881OqyzXoCwDfqtsuP\n8jw7AJeReu23Ae+ve2/r9+UDScG2irTPfq1m/ntze73553XA4fTf119P+gyszP++vmbalTTYL4Hx\n+f14FFiRl92s6Wd/kGCYkld2KqkX8OK66QcBdwKvyhvwy8Af6sOjZoc5o0lbHwPuaTDtd8C3ajb4\nSmAv0pjNeOCc/LMJMCtv9NoNeigpYMYCnyV1ScfX1PU08HZgDOk06ZoG4fFH4LD8eBKwZ5Pw+Atp\nR52Qf/92njYrv/F7AxuRgmEtzcPjWGAe8PH83LnAIfQPj22BNwMbA13AVcDxda/jFlIYvyTvQAOe\nagxhXVdS84Gn7rQlvy8Lga/m17gNcBfw1pptvhZ4V553wnqGx7uBxXXP/Qj4YX78J+ADNdOm5ffn\npa0s2yw8GmyXiaR98SOkfW83UhjParIv7wvsnH9/Nemg8a6B2svPHU7e1/P7+jhpiGBs3kce76ud\n5vvlR0kH4E1In4PXAlOa5UPT05aIWEXayQP4b6BH0i8lbZZn+RjpQ7040unIN4FdJc1stt4GpgEP\nNpj2YJ7e5xcR8fuIeI60E74HODoinoyIW0lhV/s6zoiIRyPi2Yj4HulDsX3NLPMj4qKIWAecDuzS\noI61wLaSpkVEb0Rc0+T1nBIRt0fEU6QP+675+feSjs7zI2IN6QM2lC8YnQZ8WNJU4O+Bn9e9xjsj\n4rKIeCYieoDj8ny1fhQRyyLiMdKR8pCBGhriuprZHeiKiK9HxJqIuIu0/xxcM88fI+LnEfFc3kbr\nYxLpQ1drJelUd6DpfY8nt7js+ppDCr5T8r53I6k39r6aef66L0fE0xFxZUTcnH//M+kUa6jb/kDg\njog4Pbd3NmlM5x018zTaL9eSwnXbiFgXEQvz57+hQa+25GA4PCKmAzuRuqPH58kzgRMkrZC0gtQ1\nE2lspKl8FaM3/+xDSuTNG8y+eZ7eZ1nN4y5Syi5rMB1J/yZpcb7SswLYlP5h9FDN4yeB8Q3Ow48g\npfYSSddLmtPkJdavc1J+vEVtfRHxJKl311REzCe91i8BF9Z/4CRtJukcSfdLWkXqgk6rW03tdrkn\n1/I8Q1xXMzOBLfr2i7zNvwhsVjPPsoEXHZJeUq+41hRSV3yg6X2PV7e47PqaCexRtx0+BLy8Zp76\nfXUPSVdI6pG0knSAHuq234L0vta6h/6fx0b75emkYYdzJD0g6buSxjVrbL0u1UbEElJXa6f81DLg\noxExteZnQkT8YQjr2jEiJuWfq0nn0jMkza6dT9IMYE/gt7WL1zzuAZ4Fptc8N6Nm+X2Az5MGGF8c\nEVNJRxMN6UX3r/mOiDgEeBnwHdJA28T1XM2DtbVKmkBK/KE4g3TaddoA075J2i47R8QU0qla/Wuc\nUfN4K+CBBu0MZV216ntOy4C76/aLyRHx9ibLrI9FwDaSansDu/C3QdxF9O897gI8HBGPtrjsYAba\nDr+r2w6TIuLjTZY5C/glMCMiNgV+wt+2/WDb7AFSYNXaCrh/0MIj1kbEMRExizRuMoeaQfqBDHa1\nZQdJn5U0Pf8+g9TV7euu/wQ4StKOefqmkt438NoGLf72vL4zJe0paUxe7/nAbyLiNw2WWwf8DPia\npE0k7UD/Fz2ZFC49wFhJX+X5R54hkXSopK58urQiP/3ceq7mPOAdkl4vaSPS+f9Qg+wHpLGIqwaY\nNpl01FwpaUsGvkrxCUnTJb2E1IP5nwbtDGVdtR4mjWv0uQ5YLenf8/09YyTtJGn3QdbzV5JeJGk8\nMC79qvF5e/XtKzcBR+fn300aHzg/L34acISkWfk078ukg15Lyw5B/Xa4EHilpMMkjcs/u0t6VZN1\nTAYei4in84H0gzXTekj72zYDLpkGXl8p6YOSxkr6AGmM7cLBCpf0Rkk75ys+q0inMU337cF6HquB\nPYBrJT1BCo1bSEc/IuIC0hH4nNy9vYU0sFrqk8CJpCNsL3AxaVDnPUNYblNSl+x00nli3+XkS/J6\nbid14Z6mvMt8ALBIUi/pCs7B63u+HhGLgE+RBngfJL3O5Qzh8ndEPBYRv408wlXnGOA1pF7V/5EC\ntd5ZpBH2u0gDZ8c2aGoo66p1EjArd81/ngN9Dul8+m7SKeeJpPdoqN5AutfnItLR86lce5+DSVcD\nHwe+Dbw3j88QERcD3yXdy3Ev6X0/uk3LNnMC8F5Jj0v6QUSsBt6S23uAtH9+hzTm1sg/A1+XtJo0\nHvbX2xTyKe43gN/nbb1n7YK5dzSH9Pl8lNTjnhMRtaf8jbycdGBbRbp6+jvSZ6khDbwfjm6SvgO8\nPCL+oepaBiNpEqkXs11E3D2M7SwlXQkYsAdntr5eEN9tyadXr1YymzSweUHVdTUi6R35FGsi6VLt\nzaRLqWajRkvhIWmqpPMkLclXM17XrsLW02RS1/oJ0nn894BfVFTLUBxE6sY+QLo57uAGpyJmHaul\n0xZJpwJXR8SJeTBrk4hYMdhyZjb6FYeHpE1Jo9bb+KhptuEZ6EaooXoF6dLRKZJ2Id2OfGREPFE7\nk6S5pO8/MHHixNfusMMOLTRpZs0sXLjwkYjoGom2Wul5dJMu3e4VEddKOgFYFRFfabRMd3d3LFiw\noKxSMxuUpIUR0T0SbbUyYHofcF9EXJt/P490b4CZbQCKwyMiHgKWSer7gtn+pK97m9kGoJUxD0h3\nSp6Zr7TcRfrqsZltAFoKj4i4iXSbr5ltYF4Qd5ia2chzeJhZEYeHmRVxeJhZEYeHmRVxeJhZEYeH\nmRVxeJhZEYeHmRVxeJhZEYeHmRVxeJhZkVa/VVuJRy8/f/CZRtCzd9xbdQnP8/Cfl1RdQj9n33Nb\n1SX0c33vev+HgcPqN1deUXUJ6809DzMr4vAwsyIODzMr4vAwsyIODzMr4vAwsyIODzMr4vAwsyIO\nDzMr4vAwsyIODzMr4vAwsyIODzMr4vAwsyItfSVf0lJgNbAOeDYi/P/Wmm0g2vH3PN4YEY+0YT1m\nNor4tMXMirQaHgFcKmmhpLntKMjMRodWT1v2joj7Jb0MuEzSkoi4qnaGHCpzAbbaaqsWmzOzTtFS\nzyMi7s//LgcuAGYPMM+8iOiOiO6urq5WmjOzDlIcHpImSprc9xh4C3BLuwozs87WymnLZsAFkvrW\nc1ZEXNyWqsys4xWHR0TcBezSxlrMbBTxpVozK+LwMLMiDg8zK+LwMLMiDg8zK+LwMLMiDg8zK+Lw\nMLMiDg8zK+LwMLMiDg8zK+LwMLMiDg8zK9KOP4A84pZNu63qEvqZcP9Lqy7heZ7aaOeqS+hneu8O\nVZfQz21P/qXqEkY99zzMrIjDw8yKODzMrIjDw8yKODzMrIjDw8yKODzMrIjDw8yKODzMrIjDw8yK\nODzMrIjDw8yKODzMrIjDw8yKtBweksZIulHShe0oyMxGh3b0PI4EFrdhPWY2irQUHpKmAwcCJ7an\nHDMbLVrteRwPfB54rtEMkuZKWiBpQU9PT4vNmVmnKA4PSXOA5RGxsNl8ETEvIrojorurq6u0OTPr\nMK30PPYC3ilpKXAOsJ+kM9pSlZl1vOLwiIijImJ6RGwNHAxcHhGHtq0yM+tovs/DzIq05b9eiIgr\ngSvbsS4zGx3c8zCzIg4PMyvi8DCzIg4PMyvi8DCzIg4PMyvi8DCzIg4PMyvi8DCzIg4PMyvi8DCz\nIg4PMyvSli/GjbSbJ9xSdQn9rN1un6pLeJ4Jd+9bdQn9jNeaqkvo59WP7Vx1CaOeex5mVsThYWZF\nHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFWgoP\nSeMlXSfpT5IWSTqmXYWZWWdr9e95PAPsFxG9ksYB8yX9OiKuaUNtZtbBWgqPiAigN/86Lv9Eq0WZ\nWedrecxD0hhJNwHLgcsi4tq66XMlLZC0oKenp9XmzKxDtBweEbEuInYFpgOzJe1UN31eRHRHRHdX\nV1erzZlZh2jb1ZaIWAFcARzQrnWaWedq9WpLl6Sp+fEE4M3AknYUZmadrdWrLZsDp0oaQwqicyPi\nwtbLMrNO1+rVlj8Du7WpFjMbRXyHqZkVcXiYWRGHh5kVcXiYWRGHh5kVcXiYWRGHh5kVcXiYWRGH\nh5kVcXiYWRGHh5kVcXiYWZFWv1VbiWvXPlh1Cf30brSu6hKeZ+Md11ZdQj/jHppSdQn93Ns7ueoS\nRj33PMysiMPDzIo4PMysiMPDzIo4PMysiMPDzIo4PMysiMPDzIo4PMysiMPDzIo4PMysiMPDzIo4\nPMysiMPDzIoUh4ekGZKukHSrpEWSjmxnYWbW2Vr5ex7PAp+NiBskTQYWSrosIm5tU21m1sGKex4R\n8WBE3JAfrwYWA1u2qzAz62xtGfOQtDWwG3DtANPmSlogaUFPT087mjOzDtByeEiaBJwPfCYiVtVP\nj4h5EdEdEd1dXV2tNmdmHaKl8JA0jhQcZ0bEz9pTkpmNBq1cbRFwErA4Io5rX0lmNhq00vPYCzgM\n2E/STfnn7W2qy8w6XPGl2oiYD6iNtZjZKOI7TM2siMPDzIo4PMysiMPDzIo4PMysiMPDzIo4PMys\niMPDzIo4PMysiMPDzIo4PMysiMPDzIq08jdMK3PHUyurLqGfFWMWVV3C84wZu2PVJfTTO2njqkvo\nZ+2WE6suYdRzz8PMijg8zKyIw8PMijg8zKyIw8PMijg8zKyIw8PMijg8zKyIw8PMijg8zKyIw8PM\nijg8zKyIw8PMijg8zKxIS+Eh6WRJyyXd0q6CzGx0aLXn8VPggDbUYWajTEvhERFXAY+1qRYzG0U8\n5mFmRYY9PCTNlbRA0oKenp7hbs7MRsiwh0dEzIuI7ojo7urqGu7mzGyE+LTFzIq0eqn2bOCPwPaS\n7pN0RHvKMrNO19J/vRARh7SrEDMbXXzaYmZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5m\nVsThYWZFHB5mVsThYWZFHB5mVsThYWZFWvpWbVUuee0NVZfQ+XatuoA6B1VdgLWbex5mVsThYWZF\nHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFHB5mVsThYWZFWgoP\nSQdIuk3SnZK+0K6izKzzFYeHpDHAfwFvA2YBh0ia1a7CzKyztdLzmA3cGRF3RcQa4Bz8J1/MNhit\n/CWxLYFlNb/fB+xRP5OkucDc/Oszkm5poc12mwY8UnURdTqtJtfTXKfVs/1INTTsf4YwIuYB8wAk\nLYiI7uFuc6g6rR7ovJpcT3OdWM9ItdXKacv9wIya36fn58xsA9BKeFwPbCfpFZI2Ag4Gftmessys\n0xWftkTEs5I+CVwCjAFOjohFgyw2r7S9YdJp9UDn1eR6mttg61FEjFRbZvYC4jtMzayIw8PMioxY\neHTSreySTpa0vFPuOZE0Q9IVkm6VtEjSkRXXM17SdZL+lOs5psp6+kgaI+lGSRdWXQuApKWSbpZ0\n00heIm1Sz1RJ50laImmxpNcNa3sjMeaRb2W/HXgz6Way64FDIuLWYW984HreAPQCp0XETlXUUFfP\n5sDmEXGDpMnAQuBdFW4fARMjolfSOGA+cGREXFNFPTV1/SvQDUyJiDlV1pLrWQp0R0RH3CQm6VTg\n6og4MV8B3SQiVgxXeyPV8+ioW9kj4irgsararxcRD0bEDfnxamAx6Q7equqJiOjNv47LP5WOrEua\nDhwInFhlHZ1K0qbAG4CTACJizXAGB4xceAx0K3tlH45OJmlrYDfg2orrGCPpJmA5cFlEVFoPcDzw\neeC5iuuoFcClkhbmr2FU6RVAD3BKPrU7UdLE4WzQA6YdRNIk4HzgMxGxqspaImJdROxKunN4tqTK\nTu8kzQGWR8TCqmpoYO+IeA3pm+WfyKfDVRkLvAb4cUTsBjwBDOvY4kiFh29lH0QeWzgfODMiflZ1\nPX1y1/cK4IAKy9gLeGceYzgH2E/SGRXWA0BE3J//XQ5cQDo9r8p9wH01PcTzSGEybEYqPHwrexN5\ngPIkYHFEHNcB9XRJmpofTyANdC+pqp6IOCoipkfE1qR95/KIOLSqegAkTcyD2+TTg7cAlV29i4iH\ngGWS+r5Vuz8wrAPuw/6tWii+lX3YSDob2BeYJuk+4OiIOKmqekhH1sOAm/M4A8AXI+KiiurZHDg1\nXyV7EXBuRHTE5dEOshlwQcp9xgJnRcTF1ZbEp4Az8wH6LuAjw9mYb083syIeMDWzIg4PMyvi8DCz\nIg4PMyvi8DCzIg4PMyvi8DCzIv8PkoNWkD5beMkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efe9e5920f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(111, aspect='equal')\n",
    "ax.set_xlim((0, net.shape[0]+1))\n",
    "ax.set_ylim((0, net.shape[1]+1))\n",
    "ax.set_title('Self-Organising Map after %d iterations' % n_iterations)\n",
    "\n",
    "# plot\n",
    "for x in range(1, net.shape[0] + 1):\n",
    "    for y in range(1, net.shape[1] + 1):\n",
    "        ax.add_patch(patches.Rectangle((x-0.5, y-0.5), 1, 1,\n",
    "                     facecolor=net[x-1,y-1,:],\n",
    "                     edgecolor='none'))\n",
    "plt.show()"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}