Python Programming for Biology: Bioinformatics and Beyond

Pure Python implementation

The self-organising map deals with matrices that are potentially large, and thus it would be

expected  that  a  pure  Python  implementation  will  be  pretty  slow  in  comparison  with  the

NumPy implementation. A notable reason for keeping a pure Python implementation is to

avoid the requirement of having NumPy installed, given that it is not part of the official

Python  release.  For  mainly  numerical  code  it  also  turns  out  that  a  pure  Python

implementation is fairly close in style and even in syntax to a C implementation. So once

we have a pure Python implementation it is going to be fairly easy to write the C version.

The  main  additional  burdens  in  the  C  implementation  will  be  memory  allocation,  and

wrapping Python around the C code.

The  arguments  to  the  main  function,  selfOrganisingMap,  will  be  the  same  as  for  the

NumPy  code,  except  that  we  will  use  lists  of  lists  instead  of  NumPy  arrays.  The  inputs

argument is of size the number of inputs times the depth, which in our example will be 3

because the data represents an RGB colour image. The spread argument is a list of lists of

size width times height describing how far across the self-organising map (grid) the input

influence  spreads.  The  steps  argument  is  the  number  of  iterations  for  the  algorithm.  We

need the Python random module for generating random numbers, and the math module for

taking the exponential of a number.

import math

import random

def selfOrganisingMap(inputs, spread, size, steps):

nrows, ncols = size

depth = len(inputs[0])

We then initialise the map to be random numbers between 0 and 1. We use the Python

random.random function for this. We just loop over all rows and columns and the depth.

somap = nrows * [0]

for i in range(nrows):

somap[i] = ncols * [0]

for j in range(ncols):

somap[i][j] = depth * [0]

for k in range(depth):

somap[i][j][k] = random.random()

We  then  iterate  the  number  of  times  specified,  updating  the  map  for  each  cell  (row,

column). The decay parameter determines the relative influences of the map and the input

value in updating the map.

for step in range(steps):

decay = math.exp(-step / float(steps))

---

for vector in inputs:

updateSomap(somap, vector, spread, decay)

return somap

The code that updates the map of course implements the same functionality as with the

NumPy  version  but  is  much  more  tedious  because  it  all  has  to  be  done  with  explicit

looping. First we determine the relevant sizes.

def updateSomap(somap, vector, size, spread, decay):

nrows, ncols = size

depth = len(vector)

width = len(spread)

height = len(spread[0])

Then we determine the cell with the smallest difference between the map and the input.

imin = jmin = -1

diff2min = 0

for i in range(nrows):

for j in range(ncols):

diff2 = 0.0

for k in range(depth):

diff = somap[i][j][k] - inputs[k]

diff2 += diff * diff

if (imin == -1) or (diff2 < diff2min):

imin = i

jmin = j

diff2min = diff2

Then  this  cell  and  surrounding  cells,  determined  by  the  spread  matrix,  are  updated,

using a weighting determined by the decay parameter.

halfWidth = (width-1) // 2

halfHeight = (height-1) // 2

for k in range(width):

i = (imin + k - halfWidth + nrows) % nrows

for l in range(height):

j = (jmin + l - halfHeight + ncols) % ncols

alpha = decay * spread[k][l]

for m in range(depth):

somap[i][j][m] = (1.0-alpha) * somap[i][j][m] + alpha * vector[m]