% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

addpath( '../../Code/eda_data' ); 

%close all; drawnow; 
clc; 
clear; 

% Repeating Example 3.2 - Illustrating the SMACOF algorithm for metric MDS.
% The data from this example has been saved in example32.mat.

% Use the Leukemia data, using the genes
% as the observations.
load leukemia
y = leukemia';
D = squareform(pdist(y,'seuclidean'));
[n,p] = size(D);
% Turn off this warning...
warning off MATLAB:divideByZero

% Get the first term of stress.
% This is fixed - does not depend on the configuration.
stress1 = sum(sum(D.^2))/2;  
% Now find an initial configuration - randomly generate.
% First need to compute the stress for this 2-D configuration.
d = 2;

%--
% For Exercise 3.15 we specify different random configurations: 
%--

Z = unifrnd(-2,2,n,d);
Z = unifrnd(-1,1,n,d);
Z = 2+4*randn(n,d); 

% Find the stress for this.
DZ = squareform(pdist(Z));
stress2 = sum(sum(DZ.^2))/2;
stress3 = sum(sum(D.*DZ));
oldstress = stress1 + stress2 - stress3;

% Iterate until stress converges.
tol = 10^(-6);
dstress = realmax;
numiter = 1;
dstress = oldstress;
while dstress > tol & numiter <= 100000
    numiter = numiter + 1;
    % Now get the update
    BZ = -D./DZ;
    for i = 1:n
        BZ(i,i) = 0;
        BZ(i,i) = -sum(BZ(:,i));
    end
    X = n^(-1)*BZ*Z;
    Z = X;
    % Now get the distances
    % Find the stress
    DZ = squareform(pdist(Z));
    stress2 = sum(sum(DZ.^2))/2;
    stress3 = sum(sum(D.*DZ));
    newstress = stress1 + stress2 - stress3;
    dstress = oldstress - newstress;
    oldstress = newstress;
end
    
figure
indb = strmatch('B',btcell);
indt = strmatch('T',btcell);
indna = strmatch('NA',btcell);
plot(X(indb,1),X(indb,2),'x',X(indt,1),X(indt,2),'o',X(indna,1),X(indna,2),'.')
legend({'B-cell';'T-cell';'NA'})

figure
indall = strmatch('ALL',cancertype);
indaml = strmatch('AML',cancertype);
plot(X(indall,1),X(indall,2),'x',X(indaml,1),X(indaml,2),'o')
legend({'ALL';'AML'})