%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab 5 code for making a Self Organising Feature Map grid (SOFM)
% and letting it learn the form of a letter A
%
% (c) Rajeev Raizada 
%     Dept. of Cognitive & Neural Systems
%     Boston University
%     rajeev@cns.bu.edu
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Represent the SOFM as a 2-D grid of x,y coordinates
% i.e. 3 dimensions in all: Rows, Cols, Slices
% 3rd-Dimension slice 1: x-coords 
% 3rd-Dimension slice 2: y-coords
 
num_rows = 15;
num_cols = 15;
a = 0.2;   % Exponent in eta and G reduction (Hertz, Krogh and Palmer, p.237)

%%%%% Initialise with small x and y coords, centred on the origin
dx = 0.1;
dy = 0.1;
m = cat(3,dx*(1-2*rand(num_rows,num_cols)),dy*(1-2*rand(num_rows,num_cols)));
         % The "cat" command joins up two slices along 3rd matrix dimension

%%%%% Plot the shape which will be mapped (a letter A in this case);

%%%%%%% Plot an A 
%%%%%%% NB: It's better to have this in a separate function file,
%%%%%%%     but it's harder to e.mail the program that way. 
%%%%%%%     Same for making the A input below 

figure(1);
clf;

line([-0.8 -0.4],[-1 1]);
line([-0.8 -0.5],[-1 -1]);
line([-0.5 -0.38],[-1 -0.4]);
line([-0.02  0.1],[-0.4 -1]);
line([-0.02  -0.38],[-0.4 -0.4]);
line([0 0.4],[1 -1]);
line([0 -0.4],[1 1]);
line([-0.1 -0.2],[0 0.5]);
line([-0.3 -0.2],[0 0.5]);
line([-0.3 -0.1],[0 0]);
line([0.1 0.4],[-1 -1]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The main loop

for cycle=1:5000,

	eta = cycle^(-a);		% Learning rate (how much nodes move)
	G = 0.5 + 10*cycle^(-a);		% Gaussian width parameter
	%%% NB: G<0.5 is boring because the Gaussian only covers one node

	%%%%% Give an input (in this case for an A-shape)
	x = 1-2*rand;
	y = 1-2*rand;

	while ~( (x>-0.4 & x<0 & y>-0.4 & y<0) | ... % Middle bar
	 	((y>-5*x - 0.5) & (y<-5*x + 1) ) | ...% Right diagonal
	 	((y>5*x + 1.5) & (y<5*x + 3) ) ) ...% Left diagonal
	 	x = 1-2*rand;
		y = 1-2*rand;
	end;

	inp = cat(3, x*ones(num_rows,num_cols), y*ones(num_rows,num_cols));
	      % Take the input x,y coords, and make each fill a slice
              % of a matrix the same size as m, so that they can be subtracted

	%%%%% Find winning node

	dist_mat_xy = (m - inp).^2;
              % First slice of this contains (squared) distances of x-coords,
              % second slice contains (squared) distances of y-coords
        dist_mat = sum(dist_mat_xy,3);
              % Sum across x and y slices to get total distance
	[win_rows,win_cols] = find(dist_mat==min(min(dist_mat)));
              % Finds the row and column of minimal distance grid point(s)
	rand_idx = ceil(length(win_rows)*rand);	
	win_row = win_rows(rand_idx);
	win_col = win_cols(rand_idx);
              % If two or more grid points tie for having shortest dist,
              % we need to pick one of them at random to be the winner.
              % These lines pick a random integer index, and pick the
              % entry from the winners vectors with this index
	
	%%% Calculate city-block distance from winner in grid
	[col_idx,row_idx] = meshgrid(1:num_cols,1:num_rows);
	      % This makes matrices of indices
	grid_dist = abs(row_idx-win_row) + abs(col_idx-win_col);
	
	%%% Calculate Gaussian movement-strength function for each node
	f_1dim = eta * exp(-(grid_dist/G).^2);
        f = cat(3,f_1dim,f_1dim);    % Make a slice for x and y coords

	%%%% Plot the map
	if max(cycle == [1 10 30 50 100 200 400 600 800 1000 3000 5000]),
	
		%%%%%%%%%% Do the plotting
		figure(1);
		if (cycle>1),delete(h);end;	% This wipes the old grid plot
		hold on;
		h=plot(m(:,:,1),m(:,:,2),'r-',m(:,:,1)',m(:,:,2)','r-');
				% This draws the new SOFM grid
		hold off;
		title(['Input presentation number ' num2str(cycle) ...
       			'     Neighbourhood size ' num2str(G) ...
       			'     Learning rate ' num2str(eta) ]); 
		drawnow;
		%eval(['print ' num2str(cycle) 'A.ps']);  
			% This would make a PostScript file
	end;

	%%% Move nodes
	m = m + f.*(inp-m);

end;    	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Go to next cycle