% 
% 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( '../../Code/eda_toolbox' ); 

close all; drawnow; 
clc; 
clear; 

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

% Example 3.1
% In this example, we show how to do the classical MDS. 

% First load the data - use the 'match' interpoint
% distance matrix.
load matchbpm
% First get out the data for topics 9 and 6. 
% Find the indices where they are equal to 6 or 9.
indlab = find(classlab ~= 6 & classlab ~=9);
% Now get rid of these from the distance matrix.
matchbpm(indlab,:) = [];
matchbpm(:,indlab)=[];
classlab(indlab) = [];

% Now implement the steps for classical MDS
% Find the matrix Q:
Q = -0.5*matchbpm.^2;
% Find the centering matrix H:
n = 73;
H = eye(n) - n^(-1)*ones(n,1)*ones(1,n);
% Find the matrix B:
B = H*Q*H;
% Get the eigenvalues of B.
[A,L] = eig(B);
% Find the ordering largest to smallest
[vals, inds] = sort(diag(L));
inds = flipud(inds);
vals = flipud(vals);
% Re-sort based on these.
A = A(:,inds);
L = diag(vals);

% Using 3-D for visualization purposes, 
% find the coordinates in the lower-dimensional space.
proj_dim = 3;
X = A(:,1:proj_dim)*diag(sqrt(vals(1:proj_dim)));
% Now plot in a 2-D scatterplot
ind6 = find(classlab == 6);
ind9 = find(classlab == 9);
%plot(X(ind6,1),X(ind6,2),'x',X(ind9,1),X(ind9,2),'o')
plot3(X(ind6,1),X(ind6,2),X(ind6,3),'x',X(ind9,1),X(ind9,2),X(ind9,3),'o');
xlabel('MDS Dimension 1')
ylabel('MDS Dimension 2')
zlabel('MDS Dimension 3')
legend({'Topic 6';'Topic 9'})
%saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_3_matchbpm_3d', 'epsc' );