% 
% 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' );  addpath( '../Chapter1' ); 

close all; drawnow; 
clc; 
clear; 

randn('seed',0); rand('seed',0); 

% apply some preprocessing on our data ... no 
%  
[data,midden,beachdune] = load_oronsay(0,0,0,0);

% [n=number of samples=156; p=number of features=675]
[n,p] = size(data);            

% A modified Example 3.2 - Illustrating the SMACOF algorithm for metric MDS.
%
y = data; 
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;
Z = unifrnd(-2,2,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
Xp = X; % <- pass "out" our results ... 
    
% lets plot the first two projections ... with labels ...for the beachdune labeling.
% 
ind_m = find(beachdune==0);
ind_b = find(beachdune==1); 
ind_d = find(beachdune==2); 
figure; 
h_m=plot(Xp(ind_m,1),Xp(ind_m,2),'rd','MarkerSize',8,'MarkerFaceColor','r'); hold on; 
h_b=plot(Xp(ind_b,1),Xp(ind_b,2),'gs','MarkerSize',8,'MarkerFaceColor','g');  
h_d=plot(Xp(ind_d,1),Xp(ind_d,2),'bo','MarkerSize',8,'MarkerFaceColor','b'); 
xlabel('x_1'); ylabel('x_2'); title('beachdune labeling'); 
grid on; 
legend([h_m,h_b,h_d], {'midden','beach','dune'}, 'location', 'northeast');
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_4_oronsay_MDS_BD', 'epsc' ); 

% ...for the midden labeling.
% 
ind_0 = find(midden==0);
ind_1 = find(midden==1); 
ind_2 = find(midden==2); 
figure; 
h_0=plot(Xp(ind_0,1),Xp(ind_0,2),'rd','MarkerSize',8,'MarkerFaceColor','r'); hold on; 
h_1=plot(Xp(ind_1,1),Xp(ind_1,2),'gs','MarkerSize',8,'MarkerFaceColor','g'); 
h_2=plot(Xp(ind_2,1),Xp(ind_2,2),'bo','MarkerSize',8,'MarkerFaceColor','b'); 
xlabel('x_1'); ylabel('x_2'); title('midden labeling'); 
grid on; 
legend([h_0,h_1,h_2], {'midden','Cnoc Coig','Caisteal non Gillean'}, 'location', 'northeast');
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_4_oronsay_MDS_ML', 'epsc' );