自组织映射

import numpy as np
from sklearn.datasets import load_digits
from math import sqrt
from numpy import (array,unravel_index,nditer,linalg,
random,subtract,power,exp,pi,zeros,arange,outer,meshgrid,dot)
from collections import defaultdict
from warnings import warn
from sklearn.datasets import load_digits
from pylab import plot, axis, show, pcolor, colorbar, bone
def fast_norm(x):
    return sqrt(dot(x,x.T))
class Som:
    def __init__(self,x,y,input_len,sigma=1.0,learning_rate=0.5,random_seed=None):
        if sigma>=x/2.0 or sigma>=y/2.0:
            warn('Warning: sigma is too high for the dimension of the map.')
        if random_seed:
            self.random_generator=random.RandomState(random_seed)
        else:
            self.random_generator=random.RandomState(random_seed)
        self.learning_rate=learning_rate
        self.sigma=sigma
        self.weights=self.random_generator.rand(x,y,input_len)*2-1
        self.weights=array([v/linalg.norm(v) for v in self.weights])
        self.activation_map=zeros((x,y))
        self.neigx=arange(x)
        self.neigy=arange(y)
        self.neighborhood=self.gaussian
    def _activate(self,x):
        s=subtract(x,self.weights)
        it=nditer(self.activation_map,flags=['multi_index'])
        while not it.finished:
            self.activation_map[it.multi_index]=fast_norm(s[it.multi_index])
            it.iternext()
    def activate(self,x):
        self._activate(x)
        return self.activation_map
    def gaussian(self,c,sigma):
        d=2*pi*sigma*sigma
        ax=exp(-power(self.neigx-c[0],2)/d)
        ay=exp(-power(self.neigy-c[1],2)/d)
        return outer(ax,ay)
    def diff_gaussian(self,c,sigma):
        xx,yy=meshgrid(self.neigx,self.neigy)
        p=power(xx-c[0],2)+power(yy-c[1],2)
        d=2*pi*sigma*sigma
        return exp(-(p)/d)*(1-2/d*p)
    def winner(self,x):
        self._activate(x)
        return unravel_index(self.activation_map.argmin(),self.activation_map.shape)
    def update(self,x,win,t):
        eta=self.learning_rate/(1+t/self.T)
        sig=self.sigma/(1+t/self.T)
        g=self.neighborhood(win,sig)*eta
        it=nditer(g, flags=['multi_index'])
        while not it.finished:
            self.weights[it.multi_index]+=g[it.multi_index]*(x-self.weights[it.multi_index])            
            self.weights[it.multi_index]=self.weights[it.multi_index]/fast_norm(self.weights[it.multi_index])
            it.iternext()
    def quantization(self,data):
        q=zeros(data.shape)
        for i,x in enumerate(data):
            q[i]=self.weights[self.winner(x)]
        return q
    def random_weights_init(self,data):
        it=nditer(self.activation_map, flags=['multi_index'])
        while not it.finished:
            self.weights[it.multi_index]=data[int(self.random_generator.rand()*len(data)-1)]
            self.weights[it.multi_index]=self.weights[it.multi_index]/fast_norm(self.weights[it.multi_index])
            it.iternext()
    def train_random(self,data,num_iteration):        
        self._init_T(num_iteration)        
        for iteration in range(num_iteration):
            rand_i=int(round(self.random_generator.rand()*len(data)-1))
            self.update(data[rand_i],self.winner(data[rand_i]),iteration)
    def train_batch(self,data,num_iteration):
        self._init_T(len(data)*num_iteration)
        iteration=0
        while iteration<num_iteration:
            idx=iteration%(len(data)-1)
            self.update(data[idx],self.winner(data[idx]),iteration)
            iteration+=1
    def _init_T(self,num_iteration):
        self.T=num_iteration/2 
    def distance_map(self):
        um=zeros((self.weights.shape[0],self.weights.shape[1]))
        it=nditer(um, flags=['multi_index'])
        while not it.finished:
            for ii in range(it.multi_index[0]-1,it.multi_index[0]+2):
                for jj in range(it.multi_index[1]-1,it.multi_index[1]+2):
                    if ii>=0 and ii<self.weights.shape[0] and jj>=0 and jj<self.weights.shape[1]:
                        um[it.multi_index]+=fast_norm(self.weights[ii,jj,:]-self.weights[it.multi_index])
            it.iternext()
        um=um/um.max()
        return um
    def activation_response(self,data):
        a=zeros((self.weights.shape[0],self.weights.shape[1]))
        for x in data:
            a[self.winner(x)]+=1
        return a
    def quantization_error(self,data):
        error=0
        for x in data:
            error+=fast_norm(x-self.weights[self.winner(x)])
        return error/len(data)
    def win_map(self,data):
        winmap=defaultdict(list)
        for x in data:
            winmap[self.winner(x)].append(x)
        return winmap
class TestSom:
    def setup_method(self,method):
        self.som=Som(5,5,1)
        for w in self.som.weights: 
            assert_almost_equal(1.0,linalg.norm(w))
        self.som.weights = zeros((5,5))
        self.som.weights[2,3]=5.0
        self.som.weights[1,1]=2.0
    def test_fast_norm(self):
        assert fast_norm(array([1,3]))==sqrt(1+9)
    def test_gaussian(self):
        bell=self.som.gaussian((2,2),1)
        assert bell.max()==1.0
        assert bell.argmax()==12
    def test_win_map(self):
        winners=self.som.win_map([5.0,2.0])
        assert winners[(2,3)][0]==5.0
        assert winners[(1,1)][0]==2.0
    def test_activation_reponse(self):
        response=self.som.activation_response([5.0,2.0])
        assert response[2,3]==1
        assert response[1,1]==1
    def test_activate(self):
        assert self.som.activate(5.0).argmin()==13.0 
    def test_quantization_error(self):
        self.som.quantization_error([5,2])==0.0
        self.som.quantization_error([4,1])==0.5
    def test_quantization(self):
        q=self.som.quantization(array([4,2]))
        assert q[0]==5.0
        assert q[1]==2.0
    def test_random_seed(self):
        som1=Som(5,5,2,sigma=1.0,learning_rate=0.5,random_seed=1)
        som2=Som(5,5,2,sigma=1.0,learning_rate=0.5,random_seed=1)
        assert_array_almost_equal(som1.weights,som2.weights)
        data=random.rand(100,2)
        som1=Som(5,5,2,sigma=1.0,learning_rate=0.5,random_seed=1)
        som1.train_random(data,10)
        som2=Som(5,5,2,sigma=1.0,learning_rate=0.5,random_seed=1)
        som2.train_random(data,10)
        assert_array_almost_equal(som1.weights,som2.weights)
    def test_train_batch(self):
        som=Som(5,5,2,sigma=1.0,learning_rate=0.5,random_seed=1)
        data=array([[4,2],[3,1]])
        q1=som.quantization_error(data)
        som.train_batch(data,10)
        assert q1>som.quantization_error(data)
    def test_train_random(self):
        som=Som(5,5,2,sigma=1.0,learning_rate=0.5,random_seed=1)
        data=array([[4,2],[3,1]])
        q1=som.quantization_error(data)
        som.train_random(data,10)
        assert q1>som.quantization_error(data)
    def test_random_weights_init(self):
        som=Som(2,2,2,random_seed=1)
        som.random_weights_init(array([[1.0,.0]]))
        for w in som.weights:
            assert_array_equal(w[0],array([1.0,.0]))
digits=load_digits()
data=digits.data
n_samples,n_features=data.shape
n_digits=len(np.unique(digits.target))
labels=digits.target
som=Som(16,16,64,sigma=1.3,learning_rate=0.3)
som.random_weights_init(data)
som.train_random(data,10000) 
bone()
pcolor(som.distance_map().T) 
colorbar()
labels[labels=='0']=0
labels[labels=='1']=1
labels[labels=='2']=2
labels[labels=='3']=3
labels[labels=='4']=4
labels[labels=='5']=5
labels[labels=='6']=6
labels[labels=='7']=7
labels[labels=='8']=8
labels[labels=='9']=9
markers=['o','v','1','3','8','s','p','x','D','*']
colors=["r","g","b","y","c",(0,0.1,0.8),(1,0.5,0),(1,1,0.3),"m",(0.4,0.6,0)]
for cnt,xx in enumerate(data):
 w=som.winner(xx) 
 plot(w[0]+.5,w[1]+.5,markers[labels[cnt]],markerfacecolor='None',markeredgecolor=colors[labels[cnt]],markersize=12,markeredgewidth=2)
axis([0,som.weights.shape[0],0,som.weights.shape[1]])
show()