In [1]:
'''
This code is adapted from https://codesachin.wordpress.com/2015/11/28/self-organizing-maps-with-googles-tensorflow/
'''

import tensorflow as tf
import numpy as np
 
 
class SOM(object): 
    _trained = False
 
    def __init__(self, l, m, n, dim, n_iterations=100, alpha=None, sigma=None):
        self._m = m
        self._n = n
        self._l = l
        if alpha is None:
            alpha = 0.3
        else:
            alpha = float(alpha)
        if sigma is None:
            sigma = max(l, m, n) / 2.0
        else:
            sigma = float(sigma)
        self._n_iterations = abs(int(n_iterations))
 
        self._graph = tf.Graph()
 
        with self._graph.as_default():
            self._weightage_vects = tf.Variable(tf.random_normal(
                [l*m*n, dim]))
 
            self._location_vects = tf.constant(np.array(
                list(self._neuron_locations(l, m, n))))
 
            self._vect_input = tf.placeholder("float", [dim])
            self._iter_input = tf.placeholder("float")

            bmu_index = tf.argmin(tf.sqrt(tf.reduce_sum(
                tf.pow(tf.sub(self._weightage_vects, tf.pack(
                    [self._vect_input for i in range(l*m*n)])), 2), 1)),
                                  0)
 
            slice_input = tf.pad(tf.reshape(bmu_index, [1]),
                                 np.array([[0, 1]]))
            bmu_loc = tf.reshape(tf.slice(self._location_vects, slice_input,
                                          tf.constant(np.array([1, 3]))),
                                 [3])
 
            learning_rate_op = tf.sub(1.0, tf.div(self._iter_input,
                                                  self._n_iterations))
            _alpha_op = tf.mul(alpha, learning_rate_op)
            _sigma_op = tf.mul(sigma, learning_rate_op)

            bmu_distance_squares = tf.reduce_sum(tf.pow(tf.sub(
                self._location_vects, tf.pack(
                    [bmu_loc for i in range(l*m*n)])), 2), 1)

            neighbourhood_func = tf.exp(tf.neg(tf.div(tf.cast(
                bmu_distance_squares, "float32"), tf.pow(_sigma_op, 2))))
            learning_rate_op = tf.mul(_alpha_op, neighbourhood_func)

            learning_rate_multiplier = tf.pack([tf.tile(tf.slice(
                learning_rate_op, np.array([i]), np.array([1])), [dim])
                                               for i in range(l*m*n)])
            weightage_delta = tf.mul(
                learning_rate_multiplier,
                tf.sub(tf.pack([self._vect_input for i in range(l*m*n)]),
                       self._weightage_vects))                                         
            new_weightages_op = tf.add(self._weightage_vects,
                                       weightage_delta)
            self._training_op = tf.assign(self._weightage_vects,
                                          new_weightages_op)                                       
 
            self._sess = tf.Session()
 
            init_op = tf.global_variables_initializer()
            self._sess.run(init_op)
 
    def _neuron_locations(self, l, m, n):
        for k in range(l):
            for i in range(m):
                for j in range(n):
                    yield np.array([k, i, j])
 
    def train(self, input_vects):
        for iter_no in range(self._n_iterations):
            for input_vect in input_vects:
                self._sess.run(self._training_op,
                               feed_dict={self._vect_input: input_vect,
                                          self._iter_input: iter_no})
 
        dim = len(input_vects[1])
        centroid_grid = np.zeros((self._l, self._m, self._n, dim))
        self._weightages = list(self._sess.run(self._weightage_vects))
        self._locations = list(self._sess.run(self._location_vects))
        for k in range(self._l):
            for i in range(self._m):
                for j in range(self._n):
                    centroid_grid[k,i,j,0] = self._weightages[j + (i * self._n) + (k * self._n * self._m)][0]
                    centroid_grid[k,i,j,1] = self._weightages[j + (i * self._n) + (k * self._n * self._m)][1]
                    centroid_grid[k,i,j,2] = self._weightages[j + (i * self._n) + (k * self._n * self._m)][2]
        self._centroid_grid = centroid_grid

        self._trained = True
 
    def get_centroids(self):
        if not self._trained:
            raise ValueError("SOM not trained yet")
        return self._centroid_grid
 
    def map_vects(self, input_vects): 
        if not self._trained:
            raise ValueError("SOM not trained yet")
 
        to_return = []
        for vect in input_vects:
            min_index = min([i for i in range(len(self._weightages))],
                            key=lambda x: np.linalg.norm(vect-
                                                         self._weightages[x]))
            to_return.append(self._locations[min_index])
 
        return to_return

###########################################################################
from matplotlib import pyplot as plt
 
colors = np.array(
     [[0., 0., 0.],
      [0., 0., 1.],
      [0., 0., 0.5],
      [0.125, 0.529, 1.0],
      [0.33, 0.4, 0.67],
      [0.6, 0.5, 1.0],
      [0., 1., 0.],
      [1., 0., 0.],
      [0., 1., 1.],
      [1., 0., 1.],
      [1., 1., 0.],
      [1., 1., 1.],
      [.33, .33, .33],
      [.5, .5, .5],
      [.66, .66, .66]])
color_names = \
    ['black', 'blue', 'darkblue', 'skyblue',
     'greyblue', 'lilac', 'green', 'red',
     'cyan', 'violet', 'yellow', 'white',
     'darkgrey', 'mediumgrey', 'lightgrey']
 
som = SOM(100, 20, 30, 3, 200)
som.train(colors)

image_grid = som.get_centroids()
 
mapped = som.map_vects(colors)

In [2]:
image_grid.shape


Out[2]:
(100, 20, 30, 3)

In [8]:
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()


def f(plane):
    return image_grid[plane,:,:,:]

ims = []
for i in range(100):
    im = plt.imshow(f(i), animated=True)
    ims.append([im])

ani = animation.ArtistAnimation(fig, ims, interval=100, blit=True,
                                repeat_delay=1)

ani.save('../3D_SOM_output.mp4')
#plt.show()

In [4]:
mapped


Out[4]:
[array([ 0, 19, 29]),
 array([ 9, 12, 29]),
 array([ 7, 19, 29]),
 array([ 6,  6, 26]),
 array([ 8, 12, 18]),
 array([ 9,  0, 13]),
 array([ 9, 19,  0]),
 array([0, 0, 0]),
 array([ 9,  0, 29]),
 array([9, 0, 0]),
 array([ 0, 19,  0]),
 array([ 0,  0, 21]),
 array([ 0, 19, 19]),
 array([ 1, 14, 11]),
 array([ 0,  5, 12])]

In [10]:
som.map_vects([[255,0,0]])


Out[10]:
[array([0, 0, 0])]

In [4]:
plt.imshow(image_grid[0,:,:,:])
plt.show()

Writing and loading

This next part tests writing out the neuron grid into a binary file and loading it again for redundancy.


In [80]:
with open('temp.bin', 'w') as f:
    for i in range(10):
        for j in range(20):
            for k in range(30):
                for l in range(3):
                    f.write(str(image_grid[i,j,k,l]))
                    f.write(' ')

In [83]:
with open("temp.bin", "r") as f:
    raw = f.read().split()

In [84]:
len(raw)


Out[84]:
18000

In [85]:
import numpy as np

loaded = np.zeros((10,20,30,3))
    
for i in range(10):
    for j in range(20):
        for k in range(30):
            for l in range(3):
                loaded[i,j,k,l] = float(raw[l + k*3 + j*30*3 + i*20*30*3])

In [86]:
loaded.shape


Out[86]:
(10, 20, 30, 3)

In [89]:
import matplotlib.pyplot as plt
plt.imshow(loaded[0,:,:,:])
plt.show()

In [ ]: