%MDS - Multidimensional Scaling - a variant of Sammon mapping 
%
%   [W,J,stress] = MDS(D,Y,OPTIONS)
%   [W,J,stress] = MDS(D,N,OPTIONS)
% 
% INPUT
%   D       Square (M x M) dissimilarity matrix
%   Y       M x N matrix containing starting configuration, or
%   N       Desired output dimensionality
%   OPTIONS Various parameters of the	minimization procedure put into 
%           a structure consisting of the following fields: 'q', 'optim',
%           'init','etol','maxiter', 'isratio', 'itmap', 'st' and 'inspect'
%           (default:
%             OPTIONS.q       = 0					
%             OPTIONS.optim   = 'pn'			
%             OPTIONS.init    = 'cs'
%             OPTIONS.etol    = 1e-6 (the precise value depends on q)
%             OPTIONS.maxiter = inf 
%             OPTIONS.isratio = 0 
%             OPTIONS.itmap   = 'yes' 
%             OPTIONS.st      = 1 
%             OPTIONS.inspect = 2).
% 
% OUTPUT
%   W       Multidimensional scaling mapping
%   J       Index of points removed before optimization
%   stress  Vector of stress values
%
% DESCRIPTION  
% Finds a nonlinear MDS map (a variant of the Sammon map) of objects
% represented by a symmetric distance matrix D with zero diagonal, given
% either the dimensionality N or the initial configuration Y. This is done
% in an iterative manner by minimizing the Sammon stress:
%
%   e = 1/(sum_{i<j} D_{ij}^(Q+2)) sum_{i<j} (D_{ij} - DY_{ij})^2 * D_{ij}^Q
%
% where DY is the distance matrix of Y, which should approximate D. If D(i,j) 
% = 0 for any different points i and j, then one of them is superfluous. The 
% indices of these points are returned in J.
%
% OPTIONS is an optional variable, using which the parameters for the mapping
% can be controlled. It may contain the following fields: 
%   
% 	Q				Stress measure to use (see above): -2,-1,0,1 or 2. 
%		INIT		Initialisation method for Y: 'randp', 'randnp', 'maxv', 'cs' 
%           or 'kl'. See MDS_INIT for an explanation.
% 	OPTIM		Minimization procedure to use: 'pn' for Pseudo-Newton or
%           'scg' for Scaled Conjugate Gradients.
%		ETOL		Tolerance of the minimization procedure. Usually, it should be 
%		MAXITER	in the order of 1e-6. If MAXITER is given (see below), then the 
%						optimization is stopped either when the error drops below ETOL or 
%						MAXITER iterations are reached.
% 	ISRATIO Indicates whether a ratio MDS should be performed (1) or not (0).
% 					If ISRATIO is 1, then instead of fitting the dissimilarities 
%						D_{ij}, A*D_{ij} is fitted in the stress function. The value A 
%						is estimated analytically in each iteration.
%		ITMAP		Determines the way new points are mapped, either in an iterative 
%						manner ('yes') by minimizing the stress; or by a linear projection 
%           ('no').
%		ST			Status, determines whether after each iteration the stress should 
% 	INSPECT	be printed on the screen (1) or not (0). When INSPECT > 0, 
%						ST = 1 and the mapping is onto 2D or larger, then the progress 
% 					is plotted during the minimization every INSPECT iterations.
% 
% Important:
%	1. It is assumed that D either is or approximates a Euclidean distance 
%		 matrix, i.e. D_{ij} = sqrt (sum_k(x_i - x_j)^2). 
%	2. Missing values can be handled; they should be marked by 'NaN' in D. 
%
% EXAMPLES:
% opt.optim = 'scg';
% opt.init  = 'cs'; 
% D  = sqrt(distm(a)); % Compute the Euclidean distance dataset of A
% w1 = mds(D,2,opt);   % An MDS map onto 2D initialized by Classical Scaling,
%                      % optimized by a Scaled Conjugate Gradients algorithm
% n  = size(D,1);
% y  = rand(n,2);
% w2 = mds(D,y,opt);   % An MDS map onto 2D initialized by random vectors
%
% z = rand(n,n);       % Set around 40% of the random distances to NaN, i.e. 
% z = (z+z')/2;        % not used in the MDS mapping
% z = find(z <= 0.6);
% D(z) = NaN;
% D(1:n+1:n^2) = 0;    % Set the diagonal to zero
% opt.optim = 'pn';
% opt.init  = 'randnp'; 
% opt.etol  = 1e-8;    % Should be high, as only some distances are used
% w3 = mds(D,2,opt);   % An MDS map onto 2D initialized by a random projection
%
% REFERENCES
% 1. M.F. Moler, A Scaled Conjugate Gradient Algorithm for Fast Supervised
%    Learning', Neural Networks, vol. 6, 525-533, 1993.
% 2. W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
%    Numerical Recipes in C, Cambridge University Press, Cambridge, 1992. 
% 3. I. Borg and P. Groenen, Modern Multidimensional Scaling, Springer
%    Verlag, Berlin, 1997. 
% 4. T.F. Cox and M.A.A. Cox, Multidimensional Scaling, Chapman and Hall, 
%    London, 1994.
%
% SEE ALSO
% MAPPINGS, MDS_STRESS, MDS_INIT, MDS_CS

%
% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [w,J,err,opt,y] = mds(D,y,options)

	if (nargin < 3)
		options = []; 							% Will be filled by MDS_SETOPT, below.
	end
	opt = mds_setopt(options);

	if (nargin < 2) | (isempty(y))
		prwarning(2,'no output dimensionality given, assuming N = 2.');
		y = 2; 
	end

	% No arguments given: return an untrained mapping.
	if (nargin < 1) | (isempty(D))
	  w = mapping(mfilename,{y,opt,[],[],[],[]}); 
	  w = setname(w,'MDS');
	  return
	end

	% YREP contains representative objects in the projected MDS space, i.e. 
	% for which the mapping exists. YREP is empty for the original MDS, since 
	%	no projection is available yet.

	yrep = [];              

	if (isdataset(D) | isa(D,'double'))
		[m,mm] = size(D);

		% Convert D to double, but retain labels in LAB.
		if (isdataset(D)), lab = getlab(D); D = +D; else, lab = ones(m,1); end;

		if (ismapping(y))

			% The MDS mapping exists; this means that YREP has already been stored.
			pars = getdata(y); [k,c] = size(y);

			y     = [];				% Empty, should now be found.
			yrep  = pars{1};	% There exists an MDS map, hence YREP is stored.
			opt   = pars{2};	% Options used for the mapping.
			II    = pars{3};  % Index of non-repeated points in YREP.
			winit = pars{4};	% The Classical Scaling map, if INIT = 'cs'.
			v     = pars{5};	% Weights used for adding new points if ITMAP = 'no'.
			n = c;						% Number of dimensions of the projected space.

			% Initialization by 'cs' is not possible when there is no winit 
			% (i.e. the CS map) and new points should be added.
			if (strcmp(opt.init,'cs')) & (isempty(winit))
				prwarning(2,'OPTIONS.init = cs is not possible when adding points; using kl.');
				opt.init = 'kl';	
			end

			% If YREP is a scalar, we have an empty mapping.
			if (max(size(yrep)) == 1)
				y    = yrep;
				yrep = [];
			end

			% Check whether D is a matrix with the zero diagonal for the existing map.
			if (m == mm) & (length(intersect(find(D(:)<eps),1:m+1:(m*mm))) >= m)
				w = yrep; 				% D is the same matrix as in the projection process; 
			  return						% YREP is then the solution
			end

			if (length(pars) < 6) | (isempty(pars{6}))
				yinit = [];
			else
				yinit = pars{6}; 	% Possible initial configuration for points to 
													% be added to an existing map
				if (size(yinit,1) ~= size(D,1))
					prwarning(2,'the size of the initial configuration does not match that of the dissimilarity matrix, using random initialization.')
					yinit =[];
				end
			end

		else

			% No MDS mapping available yet; perform the checks.

			if (~issym(D,1e-12))
				prwarning(2,'D is not a symmetric matrix; will average.'); 
				D = (D+D')/2;
			end

			% Check the number of zeros on the diagonal

			if (any(abs(diag(D)) > 0))
				error('D should have a zero diagonal'); 
			end
		end

	else		% D is neither a dataset nor a matrix of doubles
		error('D should be a dataset or a matrix of doubles.');
	end

	if (~isempty(y))

		% Y is the initial configuration or N, no MDS map exists yet;  
		% D is a square matrix.
														
		% Remove identical points, i.e. points for which D(i,j) = 0 for i ~= j.
		% I contains the indices of points left for the MDS mapping, J those
		% of the removed points and P those of the points left in I which were
		% identical to those removed.

		[I,J,P] = mds_reppoints(D);
		D = D(I,I); 					
		[ni,nc] = size(D);

		% NANID is an extra field in the OPTIONS structure, containing the indices
		% of NaN values (= missing values) in distance matrix D.

		opt.nanid = find(isnan(D(:)) == 1); 

		% Initialise Y.

		[m2,n] = size(y);
		if (max(m2,n) == 1)		% Y is a scalar, hence really is a dimensionality N.
			n = y;
			[y,winit] = mds_init(D,n,opt.init);
		else
			if (mm ~= m2)
				error('The matrix D and the starting configuration Y should have the same number of columns.');
			end
			winit = [];
			y = +y(I,:);
		end

		% The final number of true distances is:
		no_dist = (ni*(ni-1) - length(opt.nanid))/2;

	else											

		% This happens only if we add extra points to an existing MDS map. 
		% Remove identical points, i.e. points for which D(i,j) = 0 for i ~= j.
		% I contains the indices of points left for the MDS mapping, J those
		% of the removed points and P those of the points left in I which were
		% identical to those removed.

		[I,J,P] = mds_reppoints(D(:,II));
		D = D(I,II); 						
		[ni,nc] = size(D);
		yrep = yrep(II,:);     
		n = size(yrep,2);     

		% NANID is an extra field in the OPTIONS structure, containing the indices
		% of NaN values (= missing values) in distance matrix D.

		opt.nanid = find(isnan(D(:))); 

		% Initialise Y. if the new points should be added in an iterative manner.

		[m2,n] = size(yrep);     			
		if (~isempty(yinit)) 						% An initial configuration exists.
			y = yinit;
		elseif (strcmp(opt.init, 'cs')) & (~isempty(winit))
			y = D*winit;
		else   
			y = mds_init(D,n,opt.init);
		end

		if (~isempty(opt.nanid))       	% Rescale.
			scale = (max(yrep)-min(yrep))./(max(y)-min(y));
			y = y .* repmat(scale,ni,1); 
		end

		% The final number of true distances is:
		no_dist = (ni*nc - length(opt.nanid));
	end

	% Check whether there is enough data left.

	if (~isempty(opt.nanid))
		if (n*ni+2 > no_dist),
			error('There are too many missing distances: it is not possible to determine the MDS map.');
		end
		if (strcmp (opt.itmap,'no'))
			opt.itmap = 'yes';
			prwarning(1,'Due to the missing values, the projection can only be iterative. OPTIONS are changed appropriately.')
		end
	end

	if (opt.inspect > 0)
		opt.plotlab = lab(I,:);  % Labels to be used for plotting in MDS_SAMMAP.
	end                       

	if (isempty(yrep)) | (~isempty(yrep) & strcmp(opt.itmap, 'yes'))  

		% Either no MDS exists yet OR new points should be mapped in an iterative manner.

		printinfo(opt);
		[yy,err] = mds_sammap(D,y,yrep,opt);

		% Define the linear projection of distances.

		v = [];
		if (isempty(yrep)) & (isempty(opt.nanid))	
			if (rank(D) < m)
				v = pinv(D)*yy;
			else
				v = D \ yy;
			end
		end
	else
		% New points should be added by a linear projection of dissimilarity data.
		yy = D*v;
	end

	% Establish the projected configuration including the removed points.

	y = zeros(m,n); y(I,:) = +yy;                   
	if (~isempty(J))
		if (~isempty(yrep))
			y(J,:) = +yrep(II(P),:);
		else
			for k=length(J):-1:1				% J: indices of removed points.
				y(J(k),:) = y(P(k),:);		% P: indices of points equal to points in J.
			end 
		end
	end

	% In the definition step: shift the obtained configuration such that the
	% mean lies at the origin.

	if (isempty(yrep))
		y = y - ones(length(y),1)*mean(y);
		y = dataset(y,lab);
	else
		w = dataset(y,lab);
		return;
	end

	% In the definition step: the mapping should be stored.

	opt = rmfield(opt,'nanid');   % These fields are to be used internally only;
	opt = rmfield(opt,'plotlab'); % not to be set up from outside
	w   = mapping(mfilename,'trained',{y,opt,I,winit,v,[]},[],m,n);
	w   = setname(w,'MDS');

return

% **********************************************************************************
%                             Extra functions
% **********************************************************************************

% PRINTINFO(OPT)
%
% Prints progress information to the file handle given by OPT.ST.

function printinfo(opt)

	fprintf(opt.st,'Sammon mapping, error function with the parameter q=%d\n',opt.q);

	switch (opt.optim)
		case 'pn',  
			fprintf(opt.st,'Minimization by Pseudo-Newton algorithm\n');
		case 'scg',  
			fprintf(opt.st,'Minimization by Scaled Conjugate Gradients algorithm\n');
		otherwise 
			error(strcat('Possible initialization methods: pn (Pseudo-Newton), ',...
									 'or scg (Scaled Conjugate Gradients).'));
	end

return

% **********************************************************************************

% [I,J,P] = MDS_REPPOINTS(D)
%
% Finds the indices of repeated/left points. J contains the indices of
% repeated points in D. This means that for each j in J, there is a point
% k ~= j such that D(J(j),k) = 0. I contains the indices of the remaining 
% points, and P those of the points in I that were identical to those in J.
% Directly used in MDS routine.

function [I,J,P] = mds_reppoints(D)

	epsilon = 1e-20;			% Differences smaller than this are assumed to be zero.

	[m,mm] = size(D);

	I = 1:m; J = []; P = [];

	if (m == mm) & (all(abs(D(1:m+1:end)) <= epsilon))
		K = intersect (find (triu(ones(m),1)), find(D < epsilon));
		if (~isempty(K))
			P = mod(K,m);
			J = fix(K./m) + 1;           
			I(J) = [];                    
		end
	else
		[J,P] = find(D<=epsilon); 
		I(J)  = [];
	end

return;

% **********************************************************************************


% MDS_SETOPT  Sets parameters for the MDS mapping
%
% 	OPT = MDS_SETOPT(OPT_GIVEN)
% 
% INPUT
% 	OPT_GIVEN Parameters for the MDS mapping, described below (default: [])
%
% OUTPUT
% 	OPT 			Structure of chosen options; if OPT_GIVEN is empty, then 
% 							OPT.q       = 0
% 							OPT.optim   = 'pn'
% 							OPT.init    = 'cs'
% 							OPT.etol    = 1e-6 (the precise value depends on q)
% 							OPT.maxiter = inf 
% 							OPT.itmap   = 'yes' 
% 							OPT.isratio = 0 
% 							OPT.st      = 1 
% 							OPT.inspect = 2 
%
% DESCRIPTION  
% Parameters for the MDS mapping can be set or changed. OPT_GIVEN consists
% of the following fields: 'q', 'init', 'optim', 'etol', 'maxiter','itmap',
% 'isratio', 'st' and 'inspect'. OPTIONS can include all the fields or some
% of them only. The fields of OPT have some default values, which can be
% changed by the OPT_GIVEN field values. If OPT_GIVEN is empty, then OPT
% contains all default values. For a description of the fields, see MDS.

%
% Copyright: Elzbieta Pekalska, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function opt = mds_setopt (opt_given)

	opt.q       = 0;
	opt.init    = 'cs'; 
	opt.optim   = 'pn';
	opt.st      = 1; 
	opt.itmap   = 'yes'; 
	opt.maxiter = inf; 
	opt.inspect = 2; 
	opt.etol    = inf; 
	opt.isratio = 0; 
	opt.nanid   = [];    % Here are some extra values; set up in the MDS routine.
	opt.plotlab = [];    % Not to be changed by the user.

	if (~isempty(opt_given))
		if (~isstruct(opt_given))
			error('OPTIONS should be a structure with at least one of the following fields: q, init, etol, optim, maxiter, itmap, isratio, st or inspect.');
		end
		fn = fieldnames(opt_given);
		if (~all(ismember(fn,fieldnames(opt))))
			error('Wrong field names; valid field names are: q, init, optim, etol, maxiter, itmap, isratio, st or inspect.')
		end 
		for i = 1:length(fn)
			opt = setfield(opt,fn{i},getfield(opt_given,fn{i}));
		end
	end

	if (isempty(intersect(opt.q,-2:1:2)))
		error ('OPTIONS.q should be -2, -1, 0, 1 or 2.');
	end

	if (opt.maxiter < 2)
		error ('OPTIONS.iter should be at least 1.');
	end

	if (isinf(opt.etol))
		switch (opt.q) 								% Different defaults for different stresses.
			case -2, 
				opt.etol = 1e-6;
			case {-1,0}
				opt.etol = 10*sqrt(eps); 
			case {1,2}
				opt.etol = sqrt(eps); 
			end 
	elseif (opt.etol <= 0) | (opt.etol >= 0.1)
		error ('OPTIONS.etol should be positive and smaller than 0.1.');
	end

	if (~ismember(opt.optim, {'pn','scg'}))
		error('OPTIONS.optim should be pn or scg.');
	end

	if (~ismember(opt.itmap, {'yes','no'}))
		error('OPTIONS.itmap should be yes or no.');
	end

return



% **********************************************************************************

% MDS_SAMMAP Sammon iterative nonlinear mapping for MDS
%
% 	[YMAP,ERR] = MDS_SAMMAP(D,Y,YREP,OPTIONS)
% 
% INPUT
%   D       Square (M x M) dissimilarity matrix
%   Y       M x N matrix containing starting configuration, or
% 	YREP		Configuration of the representation points
%   OPTIONS Various parameters of the	minimization procedure put into 
%           a structure consisting of the following fields: 'q', 'optim',
%           'init','etol','maxiter', 'isratio', 'itmap', 'st' and 'inspect'
%           (default:
%             OPTIONS.q       = 0					
%             OPTIONS.optim   = 'pn'			
%             OPTIONS.init    = 'cs'
%             OPTIONS.etol    = 1e-6 (the precise value depends on q)
%             OPTIONS.maxiter = inf 
%             OPTIONS.isratio = 0 
%             OPTIONS.itmap   = 'yes' 
%             OPTIONS.st      = 1 
%             OPTIONS.inspect = 2).
% 
% OUTPUT
% 	YMAP 		Mapped configuration
% 	ERR 		Sammon stress 
%
% DESCRIPTION  
% Maps the objects given by a symmetric distance matrix D (with a zero
% diagonal) onto, say, an N-dimensional configuration YMAP by an iterative
% minimization of a variant of the Sammon stress. The minimization starts
% from the initial configuration Y; see MDS_INIT.
%
% YREP is the Sammon configuration of the representation set. It is used
% when new points have to be projected. In other words, if D is an M x M
% symmetric distance matrix, then YREP is empty; if D is an M x N matrix,
% then YMAP is sought such that D can approximate the distances between YMAP
% and YREP.
%
% Missing values can be handled by marking them by NaN in D.
%
% SEE ALSO
% MAPPINGS, MDS, MDS_CS, MDS_INIT, MDS_SETOPT

%
% Copyright: Elzbieta Pekalska, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [y,err] = mds_sammap(Ds,y,yrep,opt)

	if (nargin < 4)
		opt = []; 								% Will be filled by MDS_SETOPT, below.
	end
	if isempty(opt), opt = mds_setopt(opt); end
	if (nargin < 3)
		yrep = []; 
	end

	% Extract labels and calculate distance matrix.

	[m,n] = size(y);
	if (~isempty(yrep))
		replab = getlab(yrep);
		yrep = +yrep;
		D = sqrt(distm(y,yrep)); 
	else
		D = sqrt(distm(y));     
		yrep = [];
		replab = [];
	end

	if (isempty(opt.plotlab))
		opt.plotlab = ones(m,1);
	end

	it   = 0; 					% Iteration number.
	eold = inf; 				% Previous stress (error).

	% Calculate initial stress.

	[e,a]= mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio); err = e;

	fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e);  

	switch (opt.optim)	
		case 'pn',				% Pseudo-Newton minimization.

			epr   = e;			% Previous values of E, EOLD for line search algorithm.
			eepr  = e;
			add   = 1e-4; 	% Avoid G2 equal 0; see below.
			
			BETA 	= 1e-3; 	% Parameter for line search algorithm.
			lam1 	= 0;  	 
			lam2 	= 1;			% LAM (the step factor) lies in (LAM1, LAM2).
			LMIN 	= 1e-10;	% Minimum acceptable value for LAM in line search.
			lambest = 1e-4; % LAM = LAMBEST, if LAM was not found.

			% Loop until the error change falls below the tolerance or until
			% the maximum number of iterations is reached.

			while (abs(e-eold) >= opt.etol*(1 + abs(e))) & (it < opt.maxiter)

				% Plot progress if requested.

				if (opt.st > 0) & (opt.inspect > 0) & (mod(it,opt.inspect) == 0)
%					if (it == 0), figure(1); clf; end
					mds_plot(y,yrep,opt.plotlab,replab,e); 
				end

				yold = y; eold = e;

				% Calculate the stress.

				[e,a]      = mds_samstress (opt.q,yold,yrep,Ds,[],opt.nanid,opt.isratio);

				% Calculate gradients and pseudo-Newton update P.

				[g1,g2,cc] = mds_gradients (opt.q,yold,yrep,a*Ds,[],opt.nanid);
				p  		  = (g1/cc)./(abs(g2/cc) + add);
				slope	  = g1(:)' * p(:);

				lam = 1; 							  

				% Search for a suitable step LAM using a line search.

				while (1)

					% Take a step and calculate the delta error DE.

					y      = yold + lam .* p;  
					[e,a]  = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);
					de     = e - eold;

					% Stop if the condition for a suitable lam is fulfilled.

					if (de < BETA * lam * slope)
						break;
					end
 	
					% Try to find a suitable step LAM.

					if (lam ~= 1)
						r1 = de - lam*slope; 
						r2 = epr - eepr - lam2*slope; 
						aa = (2*de + r1/lam^2 - r2/lam2^2)/(lam2-lam1)^3;
						bb = ((-lam2*r1)/lam^2 +(lam*r2)/lam2^2)/(lam2-lam1)^2;

						if (abs(aa) <= eps)
							lamtmp = -0.5 * slope/bb;
						else
							lamtmp = (-bb + sqrt(max(0,(bb^2 - 3*aa*slope))))/(3*aa);  
						end

						% Prevent LAM from becoming too large.

						if (lamtmp > 0.5 * lam), lamtmp = 0.5 * lam; end

					else
						lamtmp = -0.5 * slope/(e - eold - slope);
					end

					% Prevent LAM from becoming too small.

					lam2 = lam;  	
					lam  = max(lamtmp, 0.1*lam);	
					epr  = e;
					eepr = eold;   

					if (lam < LMIN)
						y = yold + lambest .* p;
						[e,a] = mds_samstress(opt.q,y,yrep,Ds,[],opt.nanid,opt.isratio);
						break;
					end
				end

				it  = it + 1; 
				err = [err; e];
				fprintf(opt.st,'iteration: %4i   stress: %3.8d\n',it,e); 

			end


		case 'scg',							% Scaled Conjugate Gradient minimization.

			sigma0  = 1e-8;  										  
			lambda  = 1e-8; 			% Regularization parameter.
			lambda1 = 0;																														 

			% Calculate initial stress and direction of decrease P.

			[e,a]   = mds_samstress(opt.q,y,yrep,Ds,D,opt.nanid,opt.isratio);   
			g1      = mds_gradients(opt.q,y,yrep,a*Ds,D,opt.nanid); 	% Gradient
			p 		  = -g1;																						% Direction of decrease
			gnorm2  = g1(:)' * g1(:);
			success = 1;																							

			% Sanity check.

			if (gnorm2 < 1e-15)
				prwarning(2,['Gradient is nearly zero: ' gnorm2 ' (unlikely); initial configuration is the right one.']);  
				return
			end

			% Loop until the error change falls below the tolerance or until
			% the maximum number of iterations is reached.

			while (abs(eold - e) > opt.etol * (1 + e)) & (it < opt.maxiter)

				g0     = g1;						% Previous gradient
				pnorm2 = p(:)' * p(:);
				eold   = e;

				% Plot progress if requested.

				if (opt.inspect > 0) & (mod(it,opt.inspect) == 0) 
%					if (it == 0), figure(1); clf; end
					mds_plot(y,yrep,opt.plotlab,replab,e); 
				end

				if (success)
					sigma = sigma0/sqrt(pnorm2);			% SIGMA: a small step from y to yy
					yy    = y + sigma .*p;
					[e,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio); 
					g2    = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid); % Gradient for yy
					s     = (g2-g1)/sigma;				% Approximation of Hessian*P directly, instead of computing
					delta = p(:)' * s(:); 				% the Hessian, since only DELTA = P'*Hessian*P is in fact needed
				end

				% Regularize the Hessian indirectly to make it positive definite.
				% DELTA is now computed as P'* (Hessian + regularization * Identity) * P
				delta = delta + (lambda1 - lambda) * pnorm2;	

				% Indicate if the Hessian is negative definite; the regularization above was too small.
				if (delta < 0)
					lambda1 = 2 * (lambda - delta/pnorm2);	% Now the Hessian will be positive definite
					delta   = -delta + lambda * pnorm2; 		% This is obtained after plugging lambda1 
					lambda  = lambda1;											% into the formulation a few lines above.
				end

				mi = - p(:)' * g1(:);
				yy = y + (mi/delta) .*p;									% mi/delta is a step size 
				[ee,a] = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio); 

				% This minimization procedure is based on the second order approximation of the  
				% stress function by using the gradient and the Hessian approximation. The Hessian 
				% is regularized, but maybe not sufficiently. The ratio Dc (<=1) below indicates 
				% a proper approximation, if Dc is close to 1.

				Dc = 2 * delta/mi^2 * (e - ee); 
				e  = ee;

				% If Dc > 0, then the stress can be successfully decreased.
				success = (Dc >= 0);
				if (success)
					y       = yy;
					[ee,a]  = mds_samstress(opt.q,yy,yrep,Ds,[],opt.nanid,opt.isratio); 
					g1      = mds_gradients(opt.q,yy,yrep,a*Ds,[],opt.nanid); 
					gnorm2  = g1(:)' * g1(:);
					lambda1 = 0;

					beta = max(0,(g1(:)'*(g1(:)-g0(:)))/mi);	
					p	 = -g1 + beta .* p; 								% P is a new conjugate direction	

					if (g1(:)'*p(:) >= 0 | mod(it-1,n*m) == 0),
						p = -g1;														% No much improvement, restart
					end   

					if (Dc >= 0.75)
						lambda = 0.25 * lambda; 						% Make the regularization smaller
					end

					it  = it + 1; 
					err = [err; e];
					fprintf (opt.st,'iteration: %4i   stress: %3.8d\n',it,e); 

				else				% Dc < 0  
					% Note that for Dc < 0, the iteration number IT is not increased, 
					% so the Hessian is further regularized until SUCCESS equals 1.
					lambda1 = lambda;
				end

			% The approximation of the Hessian was poor or the stress was not 
			% decreased (Dc < 0), hence the regularization lambda is enlarged.
			if (Dc < 0.25)
				lambda = lambda + delta * (1 - Dc)/pnorm2;
			end
	  end
	end

return;



% **********************************************************************************

%MDS_STRESS - Calculate Sammon stress during optimization
%
% 	E = MDS_SAMSTRESS(Q,Y,YREP,Ds,D,NANINDEX)
%
% INPUT
% 	Q					Indicator of the Sammon stress; Q = -2,-1,0,1,2
% 	Y 				Current lower-dimensional configuration
% 	YREP			Configuration of the representation objects; it should be
% 						empty when no representation set is considered
% 	Ds				Original distance matrix
% 	D 				Approximate distance matrix (optional; otherwise computed from Y)
% 	NANINDEX	Index of the missing values; marked in Ds by NaN (optional; to
% 						  be found in Ds)
%
% OUTPUT
% 	E 				Sammon stress
%
% DESCRIPTION
% Computes the Sammon stress between the original distance matrix Ds and the
% approximated distance matrix D between the mapped configuration Y and the
% configuration of the representation set YREP, expressed as follows:
%
% E = 1/(sum_{i<j} Ds_{ij}^(q+2)) sum_{i<j} (Ds_{ij} - D_{ij})^2 * Ds_{ij}^q
%
% It is directly used in the MDS_SAMMAP routine.

%
% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [e,alpha] = mds_samstress (q,y,yrep,Ds,D,nanindex,isratio)

	% If D is not given, calculate it from Y, the current mapped points.

	if (nargin < 5) | (isempty(D))
		if (~isempty(yrep))
			D = sqrt(distm(y,yrep));     
		else
			D = sqrt(distm(y));     
		end
	end

	if (nargin < 6)
		nanindex = [];		% Not given, so calculate below.
	end

	if (nargin < 7)
		isratio = 0;			% Assume this is meant.
	end

	todefine = isempty(yrep);

	[m,n]  = size(y); [mm,k] = size(Ds);
	if (m ~= mm)
		error('The sizes of Y and Ds do not match.');
	end
	if (any(size(D) ~= size(Ds)))
		error ('The sizes of D and Ds do not match.');
	end
	m2 = m*k;

	% Convert to double.
	D  = +D; Ds = +Ds;

	% I is the index of non-NaN, non-zero (> eps) values to be included 
	% for the computation of the stress.

	I = 1:m2; 
	if (~isempty(nanindex))
		I(nanindex) = [];
	end

	O = [];
	if (todefine), O = 1:m+1:m2; end                               
  I = setdiff(I,O);

	% If OPTIONS.isratio is set, calculate optimal ALPHA to scale with.

  if (isratio)
    alpha = sum((Ds(I).^q).*D(I).^2)/sum((Ds(I).^(q+1)).*D(I));
    Ds    = alpha*Ds;
	else
		alpha = 1; 
	end
		
	% C is the normalization factor.
	c = sum(Ds(I).^(q+2)); 

	% If Q == 0, prevent unnecessary calculation (X^0 == 1).
	if (q ~= 0)
		e = sum(Ds(I).^q .*((Ds(I)-D(I)).^2))/c;
	else
		e = sum(((Ds(I)-D(I)).^2))/c;
	end

return; 



% **********************************************************************************

%MDS_GRADIENTS - Gradients for variants of the Sammon stress
%
% 	[G1,G2,CC] = MDS_GRADIENTS(Q,Y,YREP,Ds,D,NANINDEX)
%
% INPUT
% 	Q					Indicator of the Sammon stress; Q = -2,-1,0,1,2
% 	Y 				Current lower-dimensional configuration
% 	YREP			Configuration of the representation objects; it should be
% 						empty when no representation set is considered
% 	Ds				Original distance matrix
% 	D 				Approximate distance matrix (optional; otherwise computed from Y)
% 	nanindex	Index of missing values; marked in Ds by NaN (optional;
%             otherwise found from Ds)
%
% OUTPUT
% 	G1				Gradient direction
% 	G2				Approximation of the Hessian by its diagonal	
%
% DESCRIPTION  
% This is a routine used directly in the MDS_SAMMAP routine.
%
% SEE ALSO
% MDS, MDS_INIT, MDS_SAMMAP, MDS_SAMSTRESS

%
% Copyright: Elzbieta Pekalska, Robert P.W. Duin, ela@ph.tn.tudelft.nl, 2000-2003
% Faculty of Applied Sciences, Delft University of Technology
%

function [g1,g2,c] = mds_gradients(q,y,yrep,Ds,D,nanindex)

	% If D is not given, calculate it from Y, the current mapped points.

	if (nargin < 5) | (isempty(D))
		if (~isempty(yrep))
			D = sqrt(distm(y,yrep));     
		else
			D = sqrt(distm(y));     
		end
	end

	% If NANINDEX is not given, find it from Ds.

	if (nargin < 6)
		nanindex = find(isnan(Ds(:))==1);  
	end

	% YREP is empty if no representation set is defined yet.
	% This happens when the mapping should be defined.

	todefine = (isempty(yrep));
	if (todefine) 
		yrep  = y;
	end

	[m,n]  = size(y); [mm,k] = size(Ds);
	if (m ~= mm)
		error('The sizes of Y and Ds do not match.');
	end
	if (any(size(D) ~= size(Ds)))
		error ('The sizes of D and Ds do not match.');
	end
	m2 = m*k;

	% Convert to doubles.

	D  = +D; Ds = +Ds;

	% I is the index of non-NaN, non-zero (> eps) values to be included 
	% for the computation of the gradient and the Hessians diagonal.

	I = (1:m2)';
	if (~isempty(nanindex))
		I(nanindex) = [];
	end
	K    = find(Ds(I) <= eps | D(I) <= eps);
	I(K) = [];

	% C is a normalization factor.

	c = -2/sum(Ds(I).^(q+2)); 

	% Compute G1, the gradient.

	h1 = zeros(m2,1);
	if (q == 0)										% Prevent unnecessary computation when Q = 0.
		h1(I) = (Ds(I)-D(I)) ./ D(I);	
	else  
		h1(I) = (Ds(I)-D(I)) ./ (D(I).*(Ds(I).^(-q)));
	end
	h1 = reshape (h1',m,k);
	g2 = h1 * ones(k,n);					% Here G2 is assigned only temporarily,
	g1 = c * (g2.*y - h1*yrep);		% for the computation of G1.		 

	% Compute G2, the diagonal of the Hessian, if requested.

	if (nargout > 1)
		h2 = zeros(m2,1);
		switch (q)
			case -2, 
				h2(I) = -1./(Ds(I).*D(I).^3);
			case -1, 
				h2(I) = -1./(D(I).^3);
			case 0, 
				h2(I) = - Ds(I)./(D(I).^3);
			case 1, 
				h2(I) = - Ds(I).^2./(D(I).^3);
			case 2, 
				h2(I) = -(Ds(I)./D(I)).^3;
			end 
		h2 = reshape (h2',m,k);
		g2 = c * (g2 + (h2*ones(k,n)).*y.^2 + h2*yrep.^2 - 2*(h2*yrep).*y);	
	end

return;



% **********************************************************************************

%MDS_PLOT Plots the results of the MDS mapping in a 2D or 3D figure
% 
% 	MDS_PLOT(Y,YREP,LAB,REPLAB,E)
% 
% INPUT
% 	Y 			Configuration in 2D or 3D space
% 	YREP		Configuration of the representation points in 2D or 3D space
% 	LAB 		Labels of Y
% 	REPLAB	Labels of YREP
% 	E 			Stress value
%
% DESCRIPTION  
% Used directly in MDS_SAMMAP routine.


function mds_plot (y,yrep,lab,replab,e)

	% This is done in a rather complicated way in order to speed up drawing.

  K = min(size(y,2),3);
	if (K > 1)
		y   = +y;
		col = 'brmk';
		sym = ['+*xosdv^<>p']';
		ii  = [1:44];
		s   = [col(ii-floor((ii-1)/4)*4)' sym(ii-floor((ii-1)/11)*11)];
		[nlab,lablist] = renumlab([replab;lab]); c = max(nlab);
		[ms,ns] = size(s); if (ms == 1), s = setstr(ones(m,1)*s); end

		yy = [yrep; y];	

	  cla; 
	  if (K == 2)
			hold on;
			for i = 1:c
				J = find(nlab==i); plot(yy(J,1),yy(J,2),s(i,:));
			end
	  else 
			for i = 1:c
				J = find(nlab==i); plot3(yy(J,1),yy(J,2),yy(J,3),s(i,:));
				if (i == 1), hold on; grid on; end
			end
	  end

		title(sprintf('STRESS: %f', e));
	  axis equal; 
	 	drawnow;
	end

return;
