% MDS_CS Classical scaling%% 	W = MDS_CS(D,N)% % INPUT% 	D 	Square dissimilarity matrix of the size M x M% 	N 	Desired output dimensionality (optional; default: 2)%% OUTPUT% 	W		Classical scaling mapping %% DESCRIPTION  % A linear mapping of objects given by a symmetric distance matrix D with% a zero diagonal onto an N-dimensional Euclidean space such that the square % distances are preserved as much as possible. %% D is assumed to approximate the Euclidean distances, i.e. % D_{ij} = sqrt(sum_k (x_i-x_j)^2). %%% REFERENCES% 1. I. Borg and P. Groenen, Modern Multidimensional Scaling, Springer Verlag, Berlin, % 	 1997.% 2. T.F. Cox and M.A.A. Cox, Multidimensional Scaling, Chapman and Hall, London, 1994.%% SEE ALSO% MAPPINGS, 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 = mds_cs (D,n)if nargin < 2, n = 2; end	if nargin < 1 | isempty(D)		W = mapping(mfilename,n);		W = setname(W,'Classical MDS');		return	end	[m,mm] = size(D);	if (isa(D,'dataset') | isa(D,'double')) 		if isa(n,'mapping')			% The projection is already defined; apply it to new data			w = getdata(n);			[m,n] = size(D);	      if isa(D,'dataset'),				lab = getlab(D); 				D = +D;			else				lab = [];			end			D = remove_nan(D);  			D     = D.^2;			Q     = w{1}; 			me    = w{2};			L     = w{3};				Z = -repmat(1/n,n,n);                 			Z(1:n+1:end) = Z(1:n+1:end) + 1;    	 % Z = eye(n) - ones(n,n)/n			data = -0.5 * (D - me(ones(m,1),:)) * Z * Q * diag(sqrt(abs(L))./L);          			if ~isempty(lab), 				W = dataset(data,lab); 			else				W = data;			end			return; 		end	end	% Definition of the projection; we are here only if there is no projection yet	if ~issym(D,1e-12),		error('Matrix should be symmetric.')	end  D = remove_nan(D);  	D = +D.^2;	B = -repmat(1/m,m,m);						B(1:m+1:end) = B(1:m+1:end) + 1;		% B = eye(m) - ones(m,m)/m	B = -0.5 * B * D * B; 							% B is now the matrix of inner products 	ok = 0;	p  = 1;	k  = min (p*n,m);	options.disp   = 0; 	options.isreal = 1; 	options.issym  = 1; 	while ~ok & k <= m, 								% Find k largest eigenvalues		k = min (p*n,m);		[V,S,U] = eigs(B,k,'lm',options);		s = diag(real(S));		[ss,I] = sort(-s);		ok = sum(s>0) >= n;		p  = p+1;	end	if ~ok & k == m,		error ('The wanted configuration cannot be found.');	end	J = I(1:n);	W = mapping(mfilename,'trained',{V(:,J), mean(D,2)', s(J)},[],m,length(J));	W = setname(W,'Classical MDS');returnfunction D = remove_nan(D)% Remove all appearing NaN's by replacing them by the nearest neighbors.%[m,mm] = size(D);nanindex = find(isnan(D(:)));  if ~isempty(nanindex),  for i=1:m    K = find(isnan(D(i,:)));    I = 1:mm;     [md,kk] = min(D(i,:));     if md < eps,       I(kk) = [];    end    D(i,K) = min(D(i,I));  end  % check whether the diagonal is of all zeros  if length(intersect(find(D(:) < eps), 1:m+1:(m*m))) >= m,    D = (D+D')/2;  endend