|
- #!/usr/bin/python2.5
- #
- # Copyright 2009 Olivier Gillet.
- #
- # Author: Olivier Gillet (ol.gillet@gmail.com)
- #
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU General Public License as published by
- # the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- # GNU General Public License for more details.
- # You should have received a copy of the GNU General Public License
- # along with this program. If not, see <http://www.gnu.org/licenses/>.
- #
- # -----------------------------------------------------------------------------
- #
- # Self-organizing map.
-
- import numpy
- import random
-
-
- class SOM(object):
-
- def __init__(self, grid_size, radius, learning_rate):
- self._grid_size = grid_size
- self._grid_x = numpy.arange(0, grid_size * grid_size) % grid_size
- self._grid_y = numpy.arange(0, grid_size * grid_size) / grid_size
- self._radius = radius
- self._learning_rate = learning_rate
- self._codewords = None
-
- @staticmethod
- def standardize(x, std_power=1, axis=0, regularization=0.0):
- t, n = x.shape
- mean = x.mean(axis=axis)
- std = x.std(axis=axis) ** std_power + regularization
- return (x - mean) / std
-
- def classify(self, data):
- distances = ((self._codewords - data) ** 2).sum(axis=1)
- return distances.argmin(), distances
-
- def train(self, data, iterations=10000, seed=42):
- n, d = data.shape
- nodes = self._grid_size
- numpy.random.seed(seed)
- random.seed(seed)
- self._codewords = numpy.random.randn(nodes * nodes, d)
- self._error_history = []
- milestone = 2
- for i in xrange(iterations):
- if i == milestone:
- print 'iteration', i, 'of', iterations, '\t', \
- round(1000.0 * i / iterations) * 0.1, '%'
- milestone <<= 1
- radius = self._radius * 2 ** (-2 * float(i) / iterations)
- learning_rate = self._learning_rate * 2 ** (-7 * float(i) / iterations)
-
- # Pick a random vector.
- x = random.choice(data)
-
- # Find best matching unit.
- bmu, distances = self.classify(x)
- self._error_history.append(distances[bmu])
-
- # Compute neighborhood update function.
- delta_x = self._grid_x[bmu] - self._grid_x
- delta_y = self._grid_y[bmu] - self._grid_y
- rbf = numpy.exp(-(delta_x ** 2 + delta_y ** 2) / (radius * radius))
- rbf = numpy.tile(rbf.reshape((nodes * nodes, 1)), (1, d))
- update = rbf * (numpy.tile(x, (nodes * nodes, 1)) - self._codewords)
- self._codewords += learning_rate * update
- return self._error_history
-
- def checkpoint(self):
- numpy.save('weights', self._codewords)
-
- def resume(self):
- self._codewords = numpy.load('weights.npy')
-
- def plot(self, x):
- import pylab
-
- for i in xrange(x.shape[0]):
- _, d = self.classify(x[i, :])
- d = numpy.exp(-d)
- pylab.figure()
- pylab.imshow(
- d.reshape((self._grid_size, self._grid_size)),
- interpolation='nearest')
- pylab.savefig('response_%d.pdf' % i)
|
