%SAMMONM Sammon mapping
%
% W = SAMMONM(A,K,MAX)
% W = A*SAMMONM([],K,MAX)
% W = A*SAMMONM(K,MAX)
% D = B*W
%
% INPUT
% A Dataset, used for training the mapping
% B Dataset, same dimensionality as A, to be mapped
% K Target dimension of mapping (default 2)
% MAX Maximum number of iterations, default 100
%
% OUTPUT
% W Trained mapping
% D K-dimensional dataset
%
% DESCRIPTION
% This is a simplified interface to the more complex MDS routine for high
% dimensional data visualisation. The output is a non-linear projection of
% the original vector space to a K-dimensional target space.
%
% The main differences with MDS are that SAMMONM operates on feature based
% datasets, while MDS expects dissimilarity matrices; MDS maps new objects
% by a second optimization procedures minimizing the stress for the test
% objects, while SAMMONM uses a linear mapping between dissimilarities and
% the target space. See also PREX_MDS for examples. A different procedure
% for the same purpose is TSNEM.
%
% EXAMPLE
% prdatasets; % make sure prdatasets is in the path
% a = satellite; % 36D dataset, 6 classes, 6435 objects
% [x,y] = gendat(a,0.5); % split in train and test set
% w = x*sammonm; % compute mapping
% figure; scattern(x*w); % show trainset mapped to 2D: somewhat overtrained
% figure; scattern((x+randn(size(x))*1e-5)*w): % some noise helps
% figure; scattern(y*w); % show test set mapped to 2D
%
% REFERENCES
% 1. JW Sammon Jr A nonlinear mapping for data structure analysis,
% IEEE Transactions on Computers C-18, pp. 401-409,1969.
% 2. E. Pekalska, D. de Ridder, R.P.W. Duin, and M.A. Kraaijveld, A new
% method of generalizing Sammon mapping with application to algorithm
% speed-up, ASCI99, Proc. 5th Annual ASCI Conf., 1999, 221-228. [<a href="http://rduin.nl/papers/asci_99_mds.pdf">pdf</a>]
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, MAPPINGS, PCAM, MDS, TSNEM, PREX_MDS, SCATTERD, SCATTERN
% Copyright: E. Pekalska, R.P.W. Duin, r.p.w.duin@37steps.com
function out = sammonm(varargin)
argin = shiftargin(varargin,'scalar');
argin = setdefaults(argin,[],2,100,0.01);
if mapping_task(argin,'definition')
out = define_mapping(argin,'untrained');
elseif mapping_task(argin,'training')
[a,k,max_iter,pow] = deal(argin{:});
d = sqrt(distm(+a));
opt.maxiter = max_iter;
w = mds(d,a*pcam(a,k),opt);
y = getdata(w,1);
% refer to distances with a low power to get some smoothness
v = prpinv(d.^pow)*y;
out = trained_mapping(a,{v,a,pow},2);
elseif mapping_task(argin,'trained execution')
[b,w] = deal(argin{1:2});
[v,a,pow] = getdata(w);
out = setdata(b,(sqrt(distm(b,a)).^pow)*v);
else
error('Illegal call');
end
return