% 
% 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; 

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

load yeast; X = data; 
[n,p] = size(X);         % [n,p] = [number of samples=384; dimension of each sample=17]

if( 0 ) % Normalize columns (the z-scores) ... if desired 
  X = X - repmat(mean(X),size(X,1),1);
  X = X ./ repmat(std(X),size(X,1),1); 
end 

d = 2;                   % <- to compare each method we use the same target dimension for all of them 
K = 12;                  % <- the number of nearest neighbors

%
% 1) Run a non-metric multidimensional scaling algorithm
% 
fprintf('running the NMMDS algorithm...\n'); 
if( 1 ) % <- using the matlab function mdscale
  dists      = squareform(pdist(X,'cityblock'));
  [Y,stress] = mdscale(dists,d); 
else
  dists  = pdist(data,'euclidean');  % <- the dissimilarities ... 
  R = 1;                             % <- city block metric ... 
  Y = nmmds(dists,d,R);              % <- project into 10 dimensions and use the Minkowski metric (with R=1 = city block distance)
end
figure; plot(Y(:,1), Y(:,2), 'k.' );
axis tight; 
title( 'mds projection of the yeast data' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_b_yeast_mds', 'epsc' ); 

%
% 2) Now run the ISOMAP - we need the distances for input.
%
fprintf('running the ISOMAP algorithm...\n'); 
options.dims = 1:10;               % These are for ISOMAP.
options.display = 0;
[Yiso, Riso, Eiso] = isomap(dists, 'k', K, options);
% To get the coordinates, look in the structure Yiso. Note that the
% coordinates are also in the form dim x observ. If we want a 2-D
% scatter plot, then we use the 2nd record. Typing in Yiso at the
% command line shows what is in the structure.
figure; plot(Yiso.coords{2}(1,:), Yiso.coords{2}(2,:), 'k.' );
axis tight; 
title( 'isomap projection of the yeast data' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_b_yeast_isomap', 'epsc' ); 

%
% 3) Now run the LLE
%
dp = 5;                            % <- max dimension of embedding ... 
Y=lle(X.',K,dp);
figure; plot(Y(1,:),Y(2,:),'.k');
axis tight; 
title( 'lle projection of the yeast data' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_b_yeast_lle', 'epsc' ); 

%
% 4) Now run the Hessian LLE (HLLE)
%
dp = 5; 
Y = hlle(X.',K,dp);
figure; plot(Y(1,:),Y(2,:),'k.');
axis tight; 
title( 'hlle projection of the yeast data' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_b_yeast_hlle', 'epsc' ); 

%
% 5) Now run the SOM
%
fprintf('running the SOM algorithm...\n'); 
% Normalize each variable to have unit variance and mean 0.
D = som_normalize(X,'var');
sD = som_data_struct(D);
% Make the SOM
sM = som_make(sD);
sM = som_autolabel(sM,sD,'vote');
figure; 
som_show(sM,'umat','all');
som_show_add('label',sM,'subplot',1);
title( 'som projection of the yeast data' ); 
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_b_yeast_som', 'epsc' ); 

%
% 6) Now run the GTM
%
fprintf('running the GTM algorithm...\n'); 
X = data;                 % <- unnormalized data ... 
noLatPts = 400;
noBasisFn = 81;
sigma = 1.5;
[Y,MU,FI,W,beta] = gtm_stp2(X,noLatPts,noBasisFn,sigma);
lambda = 0.001;
cycles = 40;
[trndW,trndBeta,llhLog] = gtm_trn(X,FI,W,lambda,cycles,beta,'quiet');
mus = gtm_pmn(X,Y,FI,trndW,trndBeta);
%modes = gtm_pmd(X,Y,FI,trndW);
figure; plot(mus(:,1),mus(:,2),'k.');
title( 'gtm projections of the yeast data' );
saveas( gcf, '../../WriteUp/Graphics/Chapter3/prob_3_14_b_yeast_gtm', 'epsc' );