function [y, E] = sammon(x, n, opts)

%
% SAMMON - apply Sammon's nonlinear mapping 
%
%    Y = SAMMON(X) applies Sammon's nonlinear mapping procedure on
%    multivariate data X, where each row represents a pattern and each column
%    represents a feature.  On completion, Y contains the corresponding
%    co-ordinates of each point on the map.  By default, a two-dimensional
%    map is created.  Note if X contains any duplicated rows, SAMMON will
%    fail (ungracefully). 
%
%    [Y,E] = SAMMON(X) also returns the value of the cost function in E (i.e.
%    the stress of the mapping).
%
%    An N-dimensional output map is generated by Y = SAMMON(X,N) .
%
%    A set of optimisation options can also be specified using a third
%    argument, Y = SAMMON(X,N,OPTS) , where OPTS is a structure with fields:
%
%       MaxIter   - maximum number of iterations
%       TolFun    - relative tolerance on objective function
%       MaxHalves - maximum number of step halvings
%       Input     - {'raw','distance'} if set to 'distance', X is interpreted
%                   as a matrix of pairwise distances.
%       Display   - {'off', 'on', 'iter'}
%
%    The default options structure can be retrieved by calling SAMMON with
%    no parameters.
%
%    References:
%
%       [1] Sammon, John W. Jr., "A Nonlinear Mapping for Data Structure
%           Analysis", IEEE Transactions on Computers, vol. C-18, no. 5,
%           pp 401-409, May 1969.

%
% File        : sammon.m
%
% Date        : Tuesday 10th August 2004
%
% Author      : Gavin C. Cawley and Nicola L. C. Talbot
%
% Description : Simple vectorised MATLAB implementation of Sammon's non-linear
%               mapping algorithm [1].
%
% References  : [1] Sammon, John W. Jr., "A Nonlinear Mapping for Data
%                   Structure Analysis", IEEE Transactions on Computers,
%                   vol. C-18, no. 5, pp 401-409, May 1969.
%
% History     : 10/08/2004 - v1.00
%
% To do       : It might be interesting to have an option to initialise the
%               map using the first N principal components of the data (it
%               might avoid bad local minima?).
%
%               The current version does not take advantage of the symmetry
%               of the distance matrix in order to allow for easy
%               vectorisation.  This may not be a good choice for very large
%               datasets, so perhaps one day I'll get around to doing a MEX
%               version using the BLAS library etc. for very large datasets.
%
% Copyright   : (c) Dr Gavin C. Cawley, August 2004.
%
%    This program is free software; you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation; either version 2 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program; if not, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

% use the default options structure

if nargin < 3

   opts.Display   = 'iter';
   opts.Input     = 'raw';
   opts.MaxHalves = 20;
   opts.MaxIter   = 500;
   opts.TolFun    = 1e-9;

end

% the user has requested the default options structure

if nargin == 0

   y = opts;

   return;

end

% create a two-dimensional map unless dimension is specified

if nargin < 2

   n = 2;

end

% set level of verbosity

if strcmp(opts.Display, 'iter')

   display = 2;

elseif strcmp(opts.Display, 'on')

   display = 1;

else

   display = 0;

end

% create distance matrix unless given by parameters

if strcmp(opts.Input, 'distance')

   D = x;

else

   D = euclid(x, x);

end

% remaining initialisation
randn('state', 0);
N     = size(x, 1);
scale = 0.5/sum(D(:));
D     = D + eye(N);
Dinv  = 1./D;
y     = randn(N, n);
d     = euclid(y,y) + eye(N);
dinv  = 1./d;
delta = D - d;
E     = sum(sum((delta.^2).*Dinv));

% get on with it

for i=1:opts.MaxIter

   % compute gradient, Hessian and search direction

   delta = dinv - Dinv;
   g     = 4*delta*y - 4*y.*(delta*ones(N,n));
   dinv3 = dinv.^3; 
   H     = -4*delta*ones(N,n) + 4*(dinv3*(y.^2)) ...
           -8*y.*(dinv3*y) + 4*(y.^2).*(dinv3*ones(N,n));
   s     = -g(:)./abs(H(:));
   y_old = y;

   % make sure search direction is downhill

%   if s'*g(:) >= 0
%
%      s = -s; 
%
%   end

   % use step-halving procedure to ensure progress is made

   for j=1:opts.MaxHalves

      y(:) = y_old(:) + s;

      d     = euclid(y,y) + eye(N);
      dinv  = 1./d;
      delta = D - d;
      E_new = sum(sum((delta.^2).*Dinv));

      if E_new < E

         break;

      else

         s = 0.5*s;

      end

   end

   % bomb out if too many halving steps are required

   if j == opts.MaxHalves

      if display

         fprintf(1, 'Warning : MaxHalves exceeded.\n');

      end

      break;

   end

   % evaluate termination criterion

   if abs((E - E_new)/E) < opts.TolFun

      if display

         fprintf(1, 'Optimisation terminated - TolFun exceeded.\n');

      end

      break;

   end

   % report progress

   E = E_new;

   if display > 1

      fprintf(1, 'epoch = %d : E = %12.10f\n', i, E*scale);

   end

end

% fiddle stress to match the original Sammon paper

E = E*scale;

% all done

function d = euclid(x,y)

d = sqrt(sum(x.^2,2)*ones(1,size(y,1))+ones(size(x,1),1)*sum(y.^2,2)'-2*(x*y'));

