function P = som_sammon(D, varargin)

%SOM_SAMMON Computes Sammon's projection for data or map.
%
% P = som_sammon(D, odim, value, mode)
%
% ARGUMENTS ([]'s are optional) 
%
%  D        (matrix or data struct or map struct) data to be projected
%  [odim]   (scalar) output dimension, 2 by default
%  [value]  (scalar) all different modes (the next argument) require 
%            a value, default = 100
%  [mode]   (string) 'steps' or 'errlimit' or 'errchange' or 'seconds',
%            see below, default is 'steps'
%
% RETURNS   
%
%  P        (matrix) projected coordinates of each data vector
%
% The output dimension must be 2 or higher but (naturally) lower 
% than data or codebook dimension. 
%
% The mode argument determines the end condition for iteration. If 
% the mode argument is used, also the value argument has to be 
% specified. Different mode possibilities are 
%
% 'steps'      the iteration is terminated when it is run <value> 
% 'errlimit'   steps, the iteration is terminated when projection error 
%              is lower than <value>,
% 'errchange'  the iteration is terminated when change between 
%              projection error on two successive iteration rounds
%	       is less than <value> percent of total error, and
% 'seconds'    the iteration is terminated after <value> seconds 
%              of iteration.
%
% Default mode is 'steps', and default value is 100.
%
% If the first argument was a map struct, the projection matrix P 
% can be visualized using function som_showgrid.
%
% See also SOM_SHOWGRID.

% Copyright (c) 1997 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 1.0beta juuso 260997, 090997, ecco 030997, 180997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(1, 4, nargin));  % check no. of input arguments is correct

% check the type of the input argument(data matrix, data struct or map struct)

if isstruct(D)
  if isfield(D, 'data')         % data struct
    D = D.data;
  elseif isfield(D, 'codebook') % map struct
    D = D.codebook;
  else
    error('Invalid structure');
  end
end

% compute data dimensions

d     = size(D);
xsize = d(1);

if ndims(D) == 2       % D is a data matrix
  ysize = 1;
  dim   = d(2);
elseif ndims(D) == 3   % D is a codebook matrix
  ysize = d(2);
  dim   = d(3);
  D     = reshape(D, xsize * ysize, dim);
end

noc = xsize * ysize;

% output dimension

if nargin < 2 | isempty(varargin{1}), odim = 2; else odim = varargin{1}; end
if dim <= 2
  error('Will not compute Sammon''s projection for data with dimension less than 3.');
end

if odim >= dim
  error(['Projection dimension is greater than or equal to data dimension']);
elseif odim < 2
  error(['Projection dimension must be 2 or greater']);
end

% determine operating mode

if nargin > 2
  value = varargin{2};
  if nargin == 3, mode = 'steps'; else mode = varargin{3}; end
  
  if strcmp(mode, 'steps')
    runlen = value;
  elseif strcmp(mode, 'errlimit')
    errlimit = value;
  elseif strcmp(mode, 'errchange')
    errchange = value;
    e_prev = 0;
  elseif strcmp(mode, 'seconds')
    endtime = value;
  else
    error(['Illegal mode: ' mode]);
  end
else
  mode   = 'steps';     % default setting if nothing else is given
  runlen = 100;     
  fprintf(2, 'No mode was given; iterating 100 steps.\n');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialization

magic    = 0.2;          % iteration step size
noc_x_1  = ones(noc, 1); % these are used quite frequently
odim_x_1 = ones(odim,1); 

x   = rand(noc,odim);
xu  = zeros(noc,odim);

fprintf(2, 'initializing              \r');

% compute mutual distances between vectors

for i = 1:noc,
  vec                     = D(i, :);
  Tmp                     = vec(noc_x_1, 1:dim);
  Diff                    = D - Tmp; 
  Diff(find(isnan(Diff))) = 0; 
  Mdist(i, :)             = sqrt(sum((Diff.^2)'));  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% action

if strcmp(mode, 'seconds'), tic; end;
fprintf(2, 'iterating              \r');

% sammon iteration   

xd = zeros(noc, odim);
dq = zeros(noc, 1);
dr = zeros(noc, 1);

i = 0;
ready = 0;
while ~ready
  for j = 1:noc,
    xd      = -x + x(j*noc_x_1,:);
    xd2     = xd.^2;
    dpj     = sqrt(sum(xd2'))';
    dq      = Mdist(:,j) - dpj;
    dr      = Mdist(:,j) .* dpj;
    ind     = find(dr ~= 0);
    term    = dq(ind) ./ dr(ind);
    e1      = sum(xd(ind,:) .* term(:,odim_x_1));
    term2   = ((1.0 + dq(ind) ./ dpj(ind)) ./ dpj(ind)) ./ dr(ind);
    e2      = sum(term) - sum(xd2(ind,:) .* term2(:,odim_x_1));
    xu(j,:) = x(j,:) + magic * e1 ./ abs(e2);
  end

  % move the center of mass to the center 

  c = sum(xu) / noc;
  x = xu - c(noc_x_1, :);

  % compute mapping error

  e = 0; tot = 0;
  for j = 2:noc, 
    d   = Mdist(1:(j - 1), j);
    tot = tot + sum(d);
    ind = find(d ~= 0);
    xd  = -x(1:(j - 1), :) + x(j * ones(j - 1, 1), :);
    ee  = d - sqrt(sum(xd'.^2))';
    e   = e + sum(ee(ind).^2 ./ d(ind)) / tot;
  end

  i = i + 1;
  fprintf(2, '\r%d iterations, error %f', i, e);

  % determine is the iteration ready
  
  switch mode
    case 'steps', 
      if i == runlen, ready = 1; end;
    case 'errlimit',
      if e < errlimit, ready = 1; end;
    case 'errchange',
      if i > 1
	change = 100 * abs(e - e_prev) / e_prev;
	if change < errchange, ready = 1; end;
	fprintf(2, ', change of error %f %%    ', change);
      end
      e_prev = e;
    case 'seconds'
      if toc > endtime, ready = 1; end;
      fprintf(2, ', elapsed time %f seconds  ', toc);
  end
  fprintf(2, '        ');
  

  % If you want to see the Sammon's projection plotted (in 2-D and 3-D case),
  % uncomment the code below; it is not in use by default to speed up 
  % computation.
  %
  %  clf
  %  if odim == 1, 
  %    plot(x(:,1), noc_x_1, 'o');
  %  elseif odim == 2, 
  %    plot(x(:,1), x(:,2), 'o');
  %  else
  %    plot3(x(:,1), x(:,2), x(:,3), 'o')
  %  end
  %  drawnow
end

fprintf(2, '\n');

P = reshape(x, [xsize ysize odim]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%