% *********** sofm.m ********** ver. 3.2 *************** % Demonstration of Self-Organizing Feature Maps % using Kohonen's Algorithm in MATLAB % Enjoy. clc; disp(''); disp(''); disp('Formation of the Feature maps using Kohonen Algorithm ') disp(' 26 June 2000 ') disp(' Andrew P. Paplinski ') disp(' Computer Science and Software Engineering ') disp(' Monash University, Melbourne, Australia ') disp(' Copyright (c) 1997 by A.P. Paplinski ') disp(' app@csse.monash.edu.au ') disp(' ') czy = input('initialisation? Y/N [Y]: ','s'); if isempty(czy), czy = 'y' ; end if (czy == 'y') | (czy == 'Y'), clear % Generation of the input training patterns. % A two-dimensional input pattern is stored as a complex number. % First, the form of the input domain is selected: indom = menu('Select the form of the input domain:',... 'a rectangle', ... 'a triangle', ...' 'a circle', ... 'a ring' , ... 'a cross' ,... 'a letter A'); if isempty(indom), indom = 2; end % Next, the dimensionality of the output space, l, is selected. % The output units ("neurons") can be arranged in a linear, % i.e. 1-D way, or in a rectangle, i.e., in a 2-D space. el = menu('Select the dimensionality of the output domain:',... '1-dimensional output domain', ... '2-dimensional output domain'); if isempty(el), el = 1; end m1 = 12 ; m2 = 18; % m1 by m2 array of output units if (el == 1), m1 = m1*m2 ; m2 = 1 ; end m = m1*m2 ; fprintf('The output lattice is %d by %d\n', m1, m2) mOK = input('would you like to change it? Y/N [N]: ','s'); if isempty(mOK), mOK = 'n' ; end if (mOK == 'y') | (mOK == 'Y') m = 1 ; while ~((m1 > 1) & (m > 1) & (m < 4000)) m1 = input('size of the output lattice: \n m1 = ') ; if (el == 2) m2 = input('m2 = ') ; end m = m1*m2 ; end end fprintf('The output lattice is %d by %d\n', m1, m2) % The position matrix V if el == 1 V = (1:m1)' ; else [v1 v2] = meshgrid(1:m1, 1:m2); V = [v1(:) v2(:)]; end % Creating input patterns N = 20*m ; % N is the number of input vectors X = rand(1, N)+j*rand(1, N) ; ix = 1:N; if (indom == 2), ix = find((imag(X)<=2*real(X))&(imag(X)<=2-2*real(X))) ; elseif (indom == 3), ix = find(abs(X-.5*(1+j))<= 0.5) ; elseif (indom == 4), ix = find((abs(X-.5*(1+j))<= 0.5) & (abs(X-.5*(1+j)) >= 0.3)) ; elseif (indom == 5), ix = find((imag(X)<(2/3)&imag(X)>(1/3))| ... (real(X)<(2/3)&real(X)>(1/3))) ; elseif (indom == 6), ix = find((2.5*real(X)-imag(X)>0 & 2.5*real(X)-imag(X)<0.5) | ... (2.5*real(X)+imag(X)>2 & 2.5*real(X)+imag(X)<2.5) | ... (real(X)>0.2 & real(X)<0.8 & imag(X)>0.2 & imag(X)<0.4) ); end X = X(ix); N = length(X); figure(1) clf reset, hold off, % resetting workspace plot(X, '.'), title('Input Distribution') drawnow % Initialisation of weights: % a pair of real weights (one complex weight) is associated with % each output unit. Therefore the weight matrix is an m1 by m2 % matrix of complex numbers. Initial values of the weights % are taken from the input distribution : W = X(1:m).' ; X = X((m+1):N) ; N = N-m ; % as a check, the count of wins for each output unit is calculated % in the matrix "hits". The expected value of % hits = N*#epochs/m hits = zeros(m,1); % An Initial Feature Map: The map is represented by a grid spanned % by the weights of adjacent output units plotted in the input space. % In other words, the weights W(j) and W(k) are connected % by a grid segment if V(j) and V(k) are coordinates of the % adjacent output units % Initial values of the training gain, eta, and the spread, sigma, % of the neighbourhood function eta = 0.4 ; % training gain sg2i = ((m1-1)^2+(m2-1)^2)/4 ; % sg2 = 2 sigma^2 sg2 = sg2i ; figure(2) clf reset plot([0 1],[0 1],'.'), grid, hold on, if el == 1 plot(W, 'b'), else FM = full(sparse(V(:,1), V(:,2), W)) ; plot(FM, 'b'), plot(FM.', 'r') ; end title(['eta = ', num2str(eta,2), ... ' sigma^2 = ', num2str(sg2/2,3)]) hold off, drawnow % end of initialisation else % continuation eta = input('input the value of eta [0.4]: ') ; if isempty(eta), eta = 0.4; end sg2 = input(['input the value of 2sigma^2 [', ... num2str(sg2i), ']: ']) ; if isempty(sg2), sg2 = sg2i; end end % rates of diminishing the values of eta and sigma during training % We aim at final values: eta = 0.1 and sigma = 1 (sg2=2) reta = (0.2)^(2/N); rsigma = (1/sg2)^(2/N) ; % ********************** main loop ****************** frm = 1; for n = 1:N % for each input pattern, X(n), and for each output unit % which store the weight vector W(v1, v2), the distance % between X(n) and W is calculate WX = X(n) - W ; % Coordinates of the winning neuron, V(kn, :), i.e., % the neuron for which abs(WX) attains minimum [mnm kn] = min(abs(WX)); vkn = V(kn, :) ; hits(kn) = hits(kn)+1; % utilization of neurons % The neighbourhood function, NB, of the "bell" shape, % is centered around the winning unit V(kn, :) rho2 = sum(((vkn(ones(m, 1), :) - V).^2), 2) ; NB = exp(-rho2/sg2) ; % Finally, the weights are updated according to the % Kohonen learning law: W = W + eta*NB.*WX ; % Values of "eta", and "sigma" are reduced % if (n==floor(N/2)), eta = 0.8; end if (n