sammon.m

用matlab实现利用统计混沌方法解决非线性系统时间序列预测的问题

function P = sammon(D, P, varargin)

%SAMMON Computes Sammon's projection for data or map.
%
% P = sammon(D, P, [value], [mode], [alpha], [Mdist])
%
% ARGUMENTS ([]'s are optional) 
%
%  D        (matrix or data struct or map struct) data to be projected
%            NOTE: the data may _not_ have unknown components
%  P        (scalar or matrix) output dimension, or the initial
%            projection matrix, size dlen x odim
%  [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'
%  [alpha]  (scalar) iteration step size, default = 0.2
%  [Mdist]  (matrix) mutual distances between data vectors in the input 
%            space, size dlen x dlen, by default the euclidian distances
%            are used
% 
% RETURNS   
%
%  P        (matrix) projected coordinates of each data vector
%
% The output dimension must be 2 or higher but (naturally) lower 
% than data set 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.
%
% EXAMPLES
%
% P = sammon(sMap,3)  
%  projects the codebook vectors of the map to 3-dimensions. The 
%  projection can be visualized with som_showgrid like this: 
%  som_showgrid(sMap,P)
%
% P = sammon(sMap,3,[],[],[],Md)
%  same as above, but uses a precalculated mutual distance matrix Md
%
% P = sammon(sMap,som_projection(sMap,3,'pca'),20,'steps')
%  same as above, but with PCA initialization, in which case the 
%  number of training steps can be lowered
%
% P = sammon(sData,2,60,'seconds')
%  projects the data matrix of sData to a plane iterating for 60 seconds
%
%
% See also CCA, PCA, SOM_GRID.

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

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

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

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

% input data
if isstruct(D),
  if isfield(D, 'data'),         D = D.data; % data struct
  elseif isfield(D, 'codebook'), D = D.codebook; % map struct
  else error('Invalid structure');
  end
end
if any(isnan(D)), 
  error('Cannot make Sammon''s projection for data with unknown components')
end 

% compute data dimensions
orig_si = size(D); 
dim = orig_si(end); 
noc = prod(orig_si)/dim;
if length(orig_si)>2, D = reshape(D,[noc dim]); end

% output dimension / initial projection matrix
if prod(size(P))==1, 
  odim = P; 
  P = rand(noc,odim)-0.5; 
else 
  si = size(P);
  odim = si(end);
  if prod(si) ~= noc*odim, 
    error('Initial projection matrix size does not match data size');
  end
  if length(si)>2, P = reshape(P,[noc odim]); end
  inds = find(isnan(P)); 
  if length(inds), P(inds) = rand(size(inds)); end
end
if odim > dim | odim < 2, 
  error('Output dimension must be within [2, dimension of data]');
end

% determine operating mode
if nargin < 3 | isnan(varargin{1}) | isempty(varargin{1}), value=100;
else value = varargin{1}; end

if nargin < 4 | isnan(varargin{2}) | isempty(varargin{2}), mode='steps';
else mode = varargin{2}; end  
switch mode,
case 'steps',     runlen = value;
case 'errlimit',  errlimit = value;
case 'errchange', errchange = value; e_prev = 0;
case 'seconds',   endtime = value;
otherwise, error(['Illegal mode: ' mode]);
end

% iteration step size
if nargin > 4, alpha = varargin{3}; else alpha = NaN; end
if isnan(alpha) | isempty(alpha), alpha = 0.2; end

% mutual distances
if nargin > 5, Mdist = varargin{4}; else Mdist = []; end

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

% these are used quite frequently
noc_x_1  = ones(noc, 1); 
odim_x_1 = ones(odim,1); 

% compute mutual distances between vectors
if isempty(Mdist) | isnan(Mdist), 
  fprintf(2, 'computing mutual distances\r');
  for i = 1:noc,
    vec         = D(i, :);
    Diff        = D - vec(noc_x_1, 1:dim);
    Mdist(i, :) = sqrt(sum((Diff.^2)'));  
  end
else
  if size(Mdist) ~= [noc noc], 
    error('Mutual distance matrix size and data set size do not match'); 
  end
end

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

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

% sammon iteration   

x  = P ;
xu = zeros(noc, odim);
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,:) + alpha * 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
  % you can speed the computation by about 20-25% by skipping this part

  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));
  end
  e = e/tot; 

  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');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% clean up

% reshape
orig_si(end) = odim; 
P = reshape(x, orig_si);

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