wolffd@0: function net = somtrain(net, options, x)
wolffd@0: %SOMTRAIN Kohonen training algorithm for SOM.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	NET = SOMTRAIN{NET, OPTIONS, X) uses Kohonen's algorithm to train a
wolffd@0: %	SOM.  Both on-line and batch algorithms are implemented. The learning
wolffd@0: %	rate (for on-line) and neighbourhood size decay linearly. There is no
wolffd@0: %	error function minimised during training (so there is no termination
wolffd@0: %	criterion other than the number of epochs), but the  sum-of-squares
wolffd@0: %	is computed and returned in OPTIONS(8).
wolffd@0: %
wolffd@0: %	The optional parameters have the following interpretations.
wolffd@0: %
wolffd@0: %	OPTIONS(1) is set to 1 to display error values; also logs learning
wolffd@0: %	rate ALPHA and neighbourhood size NSIZE. Otherwise nothing is
wolffd@0: %	displayed.
wolffd@0: %
wolffd@0: %	OPTIONS(5) determines whether the patterns are sampled randomly with
wolffd@0: %	replacement. If it is 0 (the default), then patterns are sampled in
wolffd@0: %	order.  This is only relevant to the on-line algorithm.
wolffd@0: %
wolffd@0: %	OPTIONS(6) determines if the on-line or batch algorithm is used. If
wolffd@0: %	it is 1 then the batch algorithm is used.  If it is 0 (the default)
wolffd@0: %	then the on-line algorithm is used.
wolffd@0: %
wolffd@0: %	OPTIONS(14) is the maximum number of iterations (passes through the
wolffd@0: %	complete pattern set); default 100.
wolffd@0: %
wolffd@0: %	OPTIONS(15) is the final neighbourhood size; default value is the
wolffd@0: %	same as the initial neighbourhood size.
wolffd@0: %
wolffd@0: %	OPTIONS(16) is the final learning rate; default value is the same as
wolffd@0: %	the initial learning rate.
wolffd@0: %
wolffd@0: %	OPTIONS(17) is the initial neighbourhood size; default 0.5*maximum
wolffd@0: %	map size.
wolffd@0: %
wolffd@0: %	OPTIONS(18) is the initial learning rate; default 0.9.  This
wolffd@0: %	parameter must be positive.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	KMEANS, SOM, SOMFWD
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Check arguments for consistency
wolffd@0: errstring = consist(net, 'som', x);
wolffd@0: if ~isempty(errstring)
wolffd@0:     error(errstring);
wolffd@0: end
wolffd@0: 
wolffd@0: % Set number of iterations in convergence phase
wolffd@0: if (~options(14))
wolffd@0:     options(14) = 100;
wolffd@0: end
wolffd@0: niters = options(14);
wolffd@0: 
wolffd@0: % Learning rate must be positive
wolffd@0: if (options(18) > 0)
wolffd@0:     alpha_first = options(18);
wolffd@0: else
wolffd@0:     alpha_first = 0.9;
wolffd@0: end
wolffd@0: % Final learning rate must be no greater than initial learning rate
wolffd@0: if (options(16) > alpha_first | options(16) < 0)
wolffd@0:     alpha_last = alpha_first;
wolffd@0: else
wolffd@0:     alpha_last = options(16);
wolffd@0: end
wolffd@0: 
wolffd@0: % Neighbourhood size
wolffd@0: if (options(17) >= 0)
wolffd@0:     nsize_first = options(17);
wolffd@0: else
wolffd@0:     nsize_first = max(net.map_dim)/2;
wolffd@0: end
wolffd@0: % Final neighbourhood size must be no greater than initial size
wolffd@0: if (options(15) > nsize_first | options(15) < 0)
wolffd@0:     nsize_last = nsize_first;
wolffd@0: else
wolffd@0:     nsize_last = options(15);
wolffd@0: end
wolffd@0: 
wolffd@0: ndata = size(x, 1);
wolffd@0: 
wolffd@0: if options(6)
wolffd@0:     % Batch algorithm
wolffd@0:     H = zeros(ndata, net.num_nodes);
wolffd@0: end
wolffd@0: % Put weights into matrix form
wolffd@0: tempw = sompak(net);
wolffd@0: 
wolffd@0: % Then carry out training
wolffd@0: j = 1;
wolffd@0: while j <= niters
wolffd@0:     if options(6)
wolffd@0: 	% Batch version of algorithm
wolffd@0: 	alpha = 0.0;
wolffd@0: 	frac_done = (niters - j)/niters;
wolffd@0: 	% Compute neighbourhood
wolffd@0: 	nsize = round((nsize_first - nsize_last)*frac_done + nsize_last);
wolffd@0: 	
wolffd@0: 	% Find winning node: put weights back into net so that we can
wolffd@0: 	% call somunpak
wolffd@0: 	net = somunpak(net, tempw);
wolffd@0: 	[temp, bnode] = somfwd(net, x);
wolffd@0: 	for k = 1:ndata
wolffd@0: 	    H(k, :) = reshape(net.inode_dist(:, :, bnode(k))<=nsize, ...
wolffd@0: 		1, net.num_nodes);
wolffd@0: 	end
wolffd@0: 	s = sum(H, 1);
wolffd@0: 	for k = 1:net.num_nodes
wolffd@0: 	    if s(k) > 0
wolffd@0: 		tempw(k, :) = sum((H(:, k)*ones(1, net.nin)).*x, 1)/ ...
wolffd@0: 		    s(k);
wolffd@0: 	    end
wolffd@0: 	end
wolffd@0:     else
wolffd@0: 	% On-line version of algorithm
wolffd@0: 	if options(5)
wolffd@0: 	    % Randomise order of pattern presentation: with replacement
wolffd@0: 	    pnum = ceil(rand(ndata, 1).*ndata);
wolffd@0: 	else
wolffd@0: 	    pnum = 1:ndata;
wolffd@0: 	end
wolffd@0: 	% Cycle through dataset
wolffd@0: 	for k = 1:ndata
wolffd@0: 	    % Fraction done
wolffd@0: 	    frac_done = (((niters+1)*ndata)-(j*ndata + k))/((niters+1)*ndata);
wolffd@0: 	    % Compute learning rate
wolffd@0: 	    alpha = (alpha_first - alpha_last)*frac_done + alpha_last;
wolffd@0: 	    % Compute neighbourhood
wolffd@0: 	    nsize = round((nsize_first - nsize_last)*frac_done + nsize_last);
wolffd@0: 	    % Find best node
wolffd@0: 	    pat_diff = ones(net.num_nodes, 1)*x(pnum(k), :) - tempw;
wolffd@0: 	    [temp, bnode] = min(sum(abs(pat_diff), 2));
wolffd@0: 	
wolffd@0: 	    % Now update neighbourhood
wolffd@0: 	    neighbourhood = (net.inode_dist(:, :, bnode) <= nsize);
wolffd@0: 	    tempw = tempw + ...
wolffd@0: 		((alpha*(neighbourhood(:)))*ones(1, net.nin)).*pat_diff;
wolffd@0: 	end
wolffd@0:     end
wolffd@0:     if options(1)
wolffd@0: 	% Print iteration information
wolffd@0: 	fprintf(1, 'Iteration %d; alpha = %f, nsize = %f. ', j, alpha, ...
wolffd@0: 	nsize);
wolffd@0: 	% Print sum squared error to nearest node
wolffd@0: 	d2 = dist2(tempw, x);
wolffd@0: 	fprintf(1, 'Error = %f\n', sum(min(d2)));
wolffd@0:     end
wolffd@0:     j = j + 1;
wolffd@0: end
wolffd@0: 
wolffd@0: net = somunpak(net, tempw);
wolffd@0: options(8) = sum(min(dist2(tempw, x)));