Kohonen self organising map

import numpy as np
from matplotlib import pyplot as plt


class Kohonen(object):
    """ Here be docstring ...
    Or ideally, there would be. Sorry. There isn't. This code works, but the
    convergence is not optimised and there are loads of magic numbers. I'll
    be happy to answer questions, and a co-worker has derived a more apt
    algorithm for convergence, which I'll see if he'd be willing to share.

    Meanwhile, share and enjoy this under an MIT license!

    -- Skymandr
    skymandr sprialgalaxy fripost pinprick org
    http://eintrittverboten.wordpress.com/
    """

    def __init__(self, domain, eta_0=-1, tau_0=-1, T_o=-1, sigma_0=-1,
                 nneurons=-1, T_c=-1):
        self.domain = self.load_image(domain)

        self.tau_0 = tau_0
        self.eta_0 = eta_0
        self.T_o = T_o

        if nneurons <= 0:
            self.nneurons = self.get_nneurons()
        else:
            self.nneurons = nneurons

        if sigma_0 <= 0:
            self.sigma_0 = 0.25 * self.nneurons
        else:
            self.sigma_0 = sigma_0

        if T_o <= 0:
            self.T_o = 42 * self.nneurons
        else:
            self.T_o = T_o

        if T_c <= 0:
            self.T_c = 420 * self.nneurons
        else:
            self.T_c = T_c

        if tau_0 <= 0:
            self.tau_0 = self.T_o
        else:
            self.tau_0 = tau_0

        if eta_0 <= 0:
            ny, nx = self.domain.shape
            self.eta_0 = np.sqrt(nx * ny / self.nneurons)
        else:
            self.eta_0 = eta_0

        self.tau_1 = self.tau_0 / np.log(self.sigma_0)
        self.tau_2 = self.tau_0

        self.eta = self.eta_0
        self.sigma = self.sigma_0

        self.reset_sim()

    def get_winner(self, point, neurons):
        #Delta2 = self.get_Delta(point, neurons) ** 2).sum(0)
        DeltaChess = np.abs(self.get_Delta(point, neurons)).min(1)

        #return Delta2.argmin()
        return DeltaChess.argmin()

    def get_point(self, domain):
        thresh = np.random.random()
        ny, nx = domain.shape

        val = -1
        while val < thresh:
            x = np.random.random() * nx - 0.5
            y = np.random.random() * ny - 0.5
            // Simple probability function:
            val = domain[np.int(y), np.int(x)]
            // Better contrast:
            val = 1 - (1 - val) ** 2

        return np.array([x, y])

    def get_points(self, N):
        P = self.get_point(self.domain)

        for n in xrange(N - 1):
            P = np.vstack((P, self.get_point(self.domain)))

        return P.T

    def get_distance_penalty(self, winner, neurons, sigma):
        arg_Delta = np.arange(neurons.shape[1]) - winner

        Lambda = np.exp(-arg_Delta ** 2 / (2 * sigma ** 2)) / \
                       (np.sqrt(2 * np.pi * sigma ** 2))

        return np.array([Lambda])

    def get_Delta(self, point, neurons):
        Delta = np.array([point[0] - neurons[0], point[1] - neurons[1]])

        return Delta

    def do_step(self):
        neurons = self.neurons
        sigma   = self.sigma
        eta     = self.eta
        domain  = self.domain

        point  = self.get_point(domain)
        winner = self.get_winner(point, neurons)

        Delta = (eta * self.get_distance_penalty(winner, neurons, sigma) *
                 self.get_Delta(point, neurons))
        neurons += Delta

        if self.state == 'Ordering':
            self.t += 1
            self.sigma = self.sigma_0 * np.exp(-self.t / self.tau_1)
            self.eta   = self.eta_0 * np.exp(-self.t / self.tau_2)
            if self.t > self.T_o:
                self.t = 0
                self.state = 'Convergence'
                self.sigma = 1.0
                self.eta = 0.5 * self.eta_0
        elif self.state == 'Convergence':
            self.t += 1
            if self.t > self.T_c:
                return 42
        else:
            print "Got invalid state: {0}".format(self.state)
            raise OptionError(self.state)

        return 0

    def run_sim(self, debug=False):
        if debug:
            while self.do_step() != 42:
                if self.state == 'Ordering':
                    T = self.T_o
                elif self.state == 'Convergence':
                    T = self.T_c
                else:
                    print "Got invalid state: {0}".format(self.state)
                    raise OptionError(self.state)
                print "{0}: {1} ({2}%)".format(self.state, self.t,
                                        np.int(100 * np.float(self.t) / T))
                pass
        else:
            while self.do_step() != 42:
                pass

    def reset_sim(self):
        self.t = 0
        self.state = 'Ordering'
        self.neurons = self.get_points(self.nneurons)
        self.neurons = 0.1 * self.neurons + np.array([self.neurons.mean(1)]).T

    def load_image(self, image):
        image = plt.imread(image)

        while len(image.shape) > 2:
            image = image.mean(-1)

        image -= image.min()
        image /= image.max()
        image = (1 - image) ** 2

        return image

    def get_nneurons(self, courant=1.0):
        N = (1 - self.domain).sum()
        return np.int(np.ceil(N * courant))


class OptionError(Exception):
    r""" Error raised if invalid option is passed to a function.
    """

    def __init__(self, value):
        self.value = value

    def __str__(self):
        return repr(self.value)