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

close all; drawnow; 
clc; 
clear; 

run_mds=1;
run_isomap=1; 
run_lle=1;
run_hlle=1; 
run_som=1;
run_gtm=1;

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

[X] = load_oronsay(0,0,0,0); data = X; 
% n=number of samples=226
% p=number of features=12
[n,p] = size(X);         

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

% the initial dissimilarities ... 
dists = squareform(pdist(data,'euclidean'));  

%
% 1) Run a non-metric multidimensional scaling algorithm
% 
if( run_mds ) 
  fprintf('running the NMMDS algorithm...\n'); 
  if( 1 ) % <- using the matlab function mdscale
    [Y,stress] = mdscale(dists,d,'criterion','sstress','start','random'); 
  else
    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 oronsay data'] ); 
  saveas( gcf, ['../../WriteUp/Graphics/Chapter3/prob_3_14_f_oronsay_mds'], 'epsc' ); 
end
  
%
% 2) Now run the ISOMAP - we need the distances for input.
%
if( run_isomap ) 
  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 oronsay data'] );
  saveas( gcf, ['../../WriteUp/Graphics/Chapter3/prob_3_14_f_oronsay_isomap'], 'epsc' ); 
end 

%
% 3) Now run the LLE
%
if( run_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 oronsay data'] ); 
  saveas( gcf, ['../../WriteUp/Graphics/Chapter3/prob_3_14_f_oronsay_lle'], 'epsc' ); 
end

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

%
% 5) Now run the SOM
%
if( run_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 oronsay data'] ); 
  saveas( gcf, ['../../WriteUp/Graphics/Chapter3/prob_3_14_f_oronsay_som'], 'epsc' );  
end

%
% 6) Now run the GTM
%
if( run_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 oronsay data'] );
  saveas( gcf, ['../../WriteUp/Graphics/Chapter3/prob_3_14_f_oronsay_gtm'], 'epsc' ); 
end