%GTM  Trainable mapping fitting a Generative Topographic Mapping by EM
%
%   [W,L] = GTM (A,K,M,MAPTYPE,REG,EPS,MAXITER)
%
% INPUT
%   A       Dataset or double matrix
%   K       Vector containing number of nodes per dimension (default: [5 5], 2D map)
%   M       Vector containing number of basis functions per dimension (default: [10 10])   
%   MAPTYPE Map onto mean of posterior ('mean', default) or mode ('mode')
%   REG     Regularisation (default: 0)
%   EPS     Change in likelihood to stop training (default: 1e-5)
%   MAXITER Maximum number of iterations (default: inf)
%
% OUTPUT
%   W       GTM mapping
%   L       Likelihood
%
% DESCRIPTION
% Trains a Generative Topographic Mapping of any dimension, using the EM 
% algorithm.
%
% REFERENCES
% Bishop, C.M., Svensen, M. and Williams, C.K.I., "GTM: The Generative 
% Topographic Mapping", Neural Computation 10(1):215-234, 1998.
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% PLOTGTM, SOM, PRPLOTSOM

% (c) Dick de Ridder, 2003
% Information & Communication Theory Group
% Faculty of Electrical Engineering, Mathematics and Computer Science
% Delft University of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands

function [w,L] = gtm (a, arg2, M, mapto, reg, epsilon, maxiter)

		step = 10;

    if (nargin < 2)
        prwarning(3,'number of nodes K not given, assuming [5 5] (2D map)');
        arg2 = [5 5];
    end;
    
	if (~ismapping(arg2))
		if (nargin < 7)
		    maxiter = inf;
		end;
		if (nargin < 6)
            prwarning(3,'stopping criteria not given, assuming EPSILON = 1e-10');
            epsilon = 1e-10;
		end;
		if (nargin < 5)
            prwarning(3,'regularisation REG not given, assuming 0');
		    reg = 0;
		end;
		if (nargin < 4)
            prwarning(3,'mapping type MAPTO not given, assuming mean');
		    mapto = 'mean';
		end;
        if (nargin < 3)
            prwarning(3,'number of basis functions M not given, assuming 10');
            M = 5;
        end;
    end;
    
	if (nargin == 0) | isempty(a)
		W = prmapping(mfilename); 
		W = setname(W,'GTM');
		return; 
	end

	% Prevent annoying messages: we will tell the user about any problems.
	
	warning off;
	
		a = testdatasize(a);
    t = +a'; [d,N] = size(t); [m,k] = size(a); 

    % If we're to apply a trained mapping, unpack its parameters.
    
    if (ismapping(arg2))
        w = arg2; data = getdata(w);
        K = data{1}; M = data{2};
        W = data{3}; sigma = data{4}; mapto = data{5};
    else
        K = arg2;
    end;
    
    % Create a KK-dimensional grid with dimensions K(1), K(2), ... of 
	% D-dimensional grid points and store it as a D x KK matrix X. Do the
	% same for the basis function centers PHI_MU, with grid dimensions M(1), ...

	K = reshape(K,1,prod(size(K)));		% Turn into vectors.
	M = reshape(M,1,prod(size(M)));

	% Check: either K or M should be a vector, or both should be vectors of
	% the same length.

	if (length(K) > 1) & (length(M) == 1)
		M = M * ones(1,length(K));
	elseif (length(M) > 1) & (length(K) == 1)
		K = K * ones(1,length(M));
	elseif (length(K) ~= length(M))
		error ('Length of vectors K and M should be equal, or K and/or M should be scalar.');
	end;
	
	D = length(K); KK = prod(K); MM = prod(M);

	if (D > d)
		error ('Length of vectors K and M should be <= the data dimensionality.');
    end;
    
    x         = makegrid(K,D);				% Grid points.
	phi_mu    = makegrid(M,D);				% Basis function centers.
	phi_sigma = 2/(mean(M)-1);				% Basis function widths.

	% Pre-calculate Phi.

	for j = 1:KK
  	 	for i = 1:MM
    	    Phi(i,j) = exp(-(x(:,j)-phi_mu(:,i))'*(x(:,j)-phi_mu(:,i))/(phi_sigma^2));
	    end;
    end;
	
    if (~ismapping(arg2))  % Train mapping.

		tries = 0; retry = 1;
		while ((retry == 1) & (tries <= 5))
	
            % Initialisation.
	
		    ptx = zeros(KK,N);
			R   = (1/KK)*ones(KK,N);	
			C = prcov(t'); [U,DD] = preig(C); [dummy,ind] = sort(-diag(DD));
			W = U(:,ind(1:D))*x*prpinv(Phi);
			if (size(C,1) > D)
				sigma = sqrt(DD(D+1,D+1));
			else
				sigma = 1/mean(M);
			end;
	
            done = 0; iter = 0; retry = 0; likelihood = 1.0e20;
	
            while ((~done) & (~retry))
	
      		    iter = iter + 1;
          		done = 1;
	
		  	    factor1 = (1/(2*pi*sigma^2))^(d/2);
		  	    factor2 = (-1/(2*sigma^2));
				WPhi    = W*Phi;
	
				for i = 1:KK
					ptx(i,:) = factor1 * exp(factor2*sum((WPhi(:,i)*ones(1,N)-t).^2));
				end;
	
				if (~retry)
	
					s = sum(ptx); s(find(s==0)) = realmin;		% Prevent divide-by-zero.
					R = ptx ./ (ones(size(ptx,1),1)*s);
					G = diag(sum(R'));
	
                  	% M-step #1 (eqn. 12)

          	        W = (prinv(Phi*G*Phi' + reg*eye(MM))*Phi*R*t')';
	
                  	% M-step #1 (eqn. 13)
          	
					s = 0; WPhi = W*Phi;
					for i = 1:KK
						s = s + sum(R(i,:).*sum((WPhi(:,i)*ones(1,N)-t).^2));
					end;
                  	sigma = sqrt((1/(N*d)) * s);
	
					% Re-calculate log-likelihood.
	
                   	prev_likelihood = likelihood; likelihood = sum(log(mean(ptx)+realmin));
	
					if (rem (iter,step) == 0)
                        prwarning (10, sprintf('[%3d] L: %2.2f (change: %2.2e)', ...
                            iter, likelihood, abs ((likelihood - prev_likelihood)/likelihood)));
					end;
	
					% Continue?
	
                  	done = (abs ((likelihood - prev_likelihood)/likelihood) < epsilon);
      			    done = (done | (iter > maxiter));

      			    if (~isfinite(likelihood))
      			        prwarning(3,'Problem is poorly conditioned, retrying');
      			        retry = 1; tries = tries + 1;
      			    end;
      			    
      	        end;
			end;
        end;
	
		L = likelihood;
	
		if (~retry)
    		w = prmapping(mfilename,'trained',{K,M,W,sigma,mapto},[],k,D);
		    w = setname(w,'GTM');
	    else
            prwarning(3,'Problem is too poorly conditioned, giving up');
            prwarning(3,'Consider lowering K and M or increasing REG');
	        w = [];
	    end;
	    
    else    % Apply mapping.
    
        factor1 = (1/(2*pi*sigma^2))^(d/2);
		factor2 = (-1/(2*sigma^2));
		WPhi    = W*Phi;
	
		for i = 1:KK
		    ptx(i,:) = factor1 * exp(factor2*sum((WPhi(:,i)*ones(1,N)-t).^2));
		end;
	
		s = sum(ptx); s(find(s==0)) = realmin;		% Prevent divide-by-zero.
		R = ptx ./ (ones(size(ptx,1),1)*s);

		switch(mapto)
		    case 'mean',
		        out = (x*R)';
		    case 'mode',
		        [dummy,ind] = max(R); out = x(:,ind)';
            otherwise,
		        error ('unknown mapping type: should be mean or mode')
		 end;

		 w = setdata(a,out,getlabels(w));
    end;
    
    warning on;
    
return

% GRID = MAKEGRID (K,D)
%
% Create a KK = prod(K)-dimensional grid with dimensions K(1), K(2), ... of 
% D-dimensional uniformly spaced grid points on [0,1]^prod(K), and store it 
% as a D x KK matrix X.

function grid = makegrid (K,D)

	KK = prod(K);

	for h = 1:D
		xx{h} = 0:(1/(K(h)-1)):1;	% Support point means
	end;

	% Do that voodoo / that you do / so well...

    if (D==1)
  	    xm = xx;
    else
  	    cmd = '[';
        for h = 1:D-1, cmd = sprintf ('%sxm{%d}, ', cmd, h); end; 
        cmd = sprintf ('%sxm{%d}] = ndgrid(', cmd, D);
        for h = 1:D-1, cmd = sprintf ('%sxx{%d}, ', cmd, h); end; 
        cmd = sprintf ('%sxx{%d});', cmd, D); eval(cmd);
    end;
    
	cmd = 'mm = zeros(D, ';
	for h = 1:D-1, cmd = sprintf ('%s%d, ', cmd, K(h)); end; 
	cmd = sprintf ('%s%d);', cmd, K(D)); eval (cmd);

	for h = 1:D
		cmd = sprintf ('mm(%d,', h);
		for g = 1:D-1, cmd = sprintf ('%s:,', cmd); end; 
        cmd = sprintf ('%s:) = xm{%d};', cmd, h); eval (cmd);
	end;

	grid = reshape(mm,D,KK);

return