%Code for examples in H-plots for displaying non-metric dissimilarity matrices
%Code for Isomap and multimaps in matlab, and other manifold learning
%techniques
%h-plots, cmdscale, isoMDS, Sammon, Laub et al. in R



%%%%%%%%%%%%%%%%%%%%%
%Example 1 
%%%%%%%%%%%%%%%%%%%%

%dissimilarity matrix

d=[   0    1  5.0  7.0;  1    0 10.0 12.0;   5   10  0.0  1.5;    7   12  1.5  0.0]


%isomap
options = struct('dims',1:4,'overlay',1,'comp',1,...
         'display',1,'dijkstra',1,'verbose',1);
 [Y, R, E] = Isomap(d,'k',2,options);
 plot(Y.coords{2}(1,:), Y.coords{2}(2,:),'w.');
Y.coords{2}(1,:), Y.coords{2}(2,:)
et={'MA' 'VL' 'MO' 'SP'}
text(Y.coords{2}(1,:), Y.coords{2}(2,:),et);

% dlmwrite('aeroisomap.dat',Y.coords{2}',' ');

%multimap
%from dissimilarity to similarity
    
    s=12-d;
    s(1,1)=0;
    s(2,2)=0;
    s(3,3)=0;
    s(4,4)=0;
    
    ss = max(s ./ sum(s(:)), eps); %sum 1

 % Run the multiple-maps t-SNE model
    no_dims = 2;
    
     no_maps = 2;
    
    max_iter = 500;
    
   

    E=10000;
    %several inizialization, I take that with less error
    for i=1:50
        %mult_maps_tsne2 is the same as mult_maps_tsne but return the error
        %in C
[mapsi, weightsi,C] = mult_maps_tsne2(ss, no_maps, no_dims, max_iter);
if C<E
    maps=mapsi;
    weights=weightsi;
    end
    end
    
    wi=0.00000001; %small importance weight, this can be changed, it is too small
    
    for m=1:size(maps, 3)
        
        % Remove points with small importance weight
        figure(m);
        if size(maps, 3) > 1
            ii = find(weights(:,m) > wi);
        else
            ii = 1:size(maps, 1);
        end
        
        % Make scatterplot
        scatter(maps(ii, 1, m), maps(ii, 2, m), weights(ii, m) * 40,'k'); axis off
        XLim = get(gca, 'XLim');
        width = XLim(2) - XLim(1);
            
        % Place the text labels
        for i=1:length(ii)
            text(maps(ii(i), 1, m) + 0.006 * width, maps(ii(i), 2, m), et{ii(i)});
        end
        axis off
        set(gcf, 'Position', [1 + 3 * m -150 + 3 * m 1280 705]);
        drawnow
    end    
    
maps, weights

%dlmwrite('aeromulti1.dat',maps(:,:,1),' ');
%dlmwrite('aeromulti2.dat',maps(:,:,2),' ');
%dlmwrite('aeromultiw.dat',weights,' ');


%%%%%%%%%%%%%%%%%%%%%
%Example 2: : suppl. material 
%%%%%%%%%%%%%%%%%%%%

%multimap
 d1=[  0.00 2.36 2.59 1.78 4.74 4.82 4.98 4.72; 2.36 0.00 2.39 1.60 4.98 5.06 5.22 4.96; ...
 2.59 2.39 0.00 2.09 5.29 5.37 5.53 5.27; 1.78 1.60 2.09 0.00 5.08 5.16 5.32 5.06; ...
 4.74 4.98 5.29 5.08 0.00 1.20 1.82 1.62; 4.82 5.06 5.37 5.16 1.20 0.00 2.98 1.78; ...
4.98 5.22 5.53 5.32 1.82 2.98 0.00 2.02; 4.72 4.96 5.27 5.06 1.62 1.78 2.02 0.00];
 d2=[ 0.00 4.15 2.03 4.14 1.26 4.33 0.69 4.85; 4.15 0.00 4.70 0.57 4.37 1.82 4.24 2.02; ...
 2.03 4.70 0.00 4.69 1.85 4.88 1.68 5.40;  4.14 0.57 4.69 0.00 4.36 1.83 4.23 2.67; ...
 1.26 4.37 1.85 4.36 0.00 4.55 0.73 5.07; 4.33 1.82 4.88 1.83 4.55 0.00 4.42 2.14; ...
 0.69 4.24 1.68 4.23 0.73 4.42 0.00 4.94;  4.85 2.02 5.40 2.67 5.07 2.14 4.94 0.00];

d=d1-d2;
%max(d) es 4.29

s=4.29-d; %similarity

for i=1:8
    s(i,i)=0;
end


     ss = max(s ./ sum(s(:)), eps);

 % Run the multiple-maps t-SNE model
    no_dims = 2;
    
        no_maps = 2;
    
         
          max_iter = 500;
    E=10000;
    %several inizialization, I take that with less error
    for i=1:50
        %mult_maps_tsne2 is the same as mult_maps_tsne but return the error
        %in C
[mapsi, weightsi,C] = mult_maps_tsne2(ss, no_maps, no_dims, max_iter);
if C<E
    maps=mapsi;
    weights=weightsi;
    end
    end
    
     et={'1','2','3', '4', '5', '6', '7', '8'}
     
     
        wi=0.001; %small importance weight, this can be changed, it is too small
    
    for m=1:size(maps, 3)
        
        % Remove points with small importance weight
        figure(m);
        if size(maps, 3) > 1
            ii = find(weights(:,m) > wi);
        else
            ii = 1:size(maps, 1);
        end
        
        % Make scatterplot
        scatter(maps(ii, 1, m), maps(ii, 2, m), weights(ii, m) * 40); axis off
        XLim = get(gca, 'XLim');
        width = XLim(2) - XLim(1);
            
        % Place the text labels
        for i=1:length(ii)
            text(maps(ii(i), 1, m) + 0.006 * width, maps(ii(i), 2, m), et{ii(i)});
        end
        axis off
        set(gcf, 'Position', [1 + 3 * m -150 + 3 * m 1280 705]);
        drawnow
    end    
    
maps, weights

     
    %dlmwrite('laubmulti1.dat',maps(:,:,1),' ');
    %dlmwrite('laubmulti2.dat',maps(:,:,2),' ');
    %dlmwrite('laubmultiw.dat',weights,' ');

    
%%%%%%%%%%%%%%%%%%%%%
%Example 1 with the MANIfold Learning techniques
%%%%%%%%%%%%%%%%%%%%    

%dissimilarity matrix

d=[   0    1  5.0  7.0;  1    0 10.0 12.0;   5   10  0.0  1.5;    7   12  1.5  0.0];

X=d;
K=2; %nearest neighbors
di=2; %2-dim

% --- LLE ALGORITHM (using K nearest neighbors)
% Written by Sam Roweis & Lawrence Saul

[Y] = lle(d,K,di);

plot(Y(1,:), Y(2,:),'k.');
set(gca,'FontSize',28);

et={'MA' 'VL' 'MO' 'SP'}
text(Y(1,:), Y(2,:),et,'FontSize',28);

% dlmwrite('aerolle.dat',Y',' ');

%% --- HLLE function
% Written by David Donoho & Carrie Grimes, 2003.
 [Y, mse] = HLLE(X,K,di); %error, not computed
 
 %% --- leigs function for Laplacian eigenmap.
% Written by Belkin & Niyogi, 2002.
 
[Y,V] = leigs(X, 'nn',K, di+1);

plot(Y(:,1), Y(:,2),'k.');
set(gca,'FontSize',28);
et={'MA' 'VL' 'MO' 'SP'}
text(Y(:,1), Y(:,2),et,'FontSize',28);

%% % --- diffusionKernel function
% Written by R. Coifman & S. Lafon.
%[Y] = diffusionKernel (X,sigmaK,alpha,di)
[Y] = diffusionKernel (X,10,1,di);

plot(Y(:,1), Y(:,2),'k.');
set(gca,'FontSize',28);
et={'MA' 'VL' 'MO' 'SP'}
text(Y(:,1), Y(:,2),et,'FontSize',28);

% --- LTSA function
% Written by Zhenyue Zhang & Hongyuan Zha, 2004.
% Reference: http://epubs.siam.org/sam-bin/dbq/article/41915
 [Y,NI] = ltsa(X,di,K);

 plot(Y(1,:), Y(2,:),'k.');
set(gca,'FontSize',28);
et={'MA' 'VL' 'MO' 'SP'}
text(Y(1,:), Y(2,:),et,'FontSize',28);


%%%%%%%%%%%%%%%%%%%%%
%Example 3 with a = 2: suppl. material
%%%%%%%%%%%%%%%%%%%%

%isomap
n=5;
a=2;
 d=ones(5,5);

d(1,1)=0;
d(2,2)=0;
d(3,3)=0;
d(4,4)=0;
d(5,5)=0;
d(1,2)=a;
d(2,1)=a;
d(5,4)=4 ;
d(4,5)=d(5,4);
d(5,3)=2;
d(3,5)=d(5,3);

options = struct('dims',1:4,'overlay',1,'comp',1,...
         'display',1,'dijkstra',1,'verbose',1);
 [Y, R, E] = Isomap(d,'k',2,options);
 plot(Y.coords{2}(1,:), Y.coords{2}(2,:),'w.');
Y.coords{2}(1,:), Y.coords{2}(2,:)
et={'A' 'B' 'C' 'D' 'E'};
text(Y.coords{2}(1,:), Y.coords{2}(2,:),et);

% dlmwrite('brand2isomap.dat',Y.coords{2}',' ');



%%%%%%%%%%%%%%%%%%%%%
%Example 3 with the MANIfold Learning techniques
%%%%%%%%%%%%%%%%%%%%    

n=5;
a=2;
 d=ones(5,5);

d(1,1)=0;
d(2,2)=0;
d(3,3)=0;
d(4,4)=0;
d(5,5)=0;
d(1,2)=a;
d(2,1)=a;
d(5,4)=4 ;
d(4,5)=d(5,4);
d(5,3)=2;
d(3,5)=d(5,3);


et={'A' 'B' 'C' 'D' 'E'};

X=d;
K=2; %nearest neighbors
di=2; %2-dim

% --- LLE ALGORITHM (using K nearest neighbors)
% Written by Sam Roweis & Lawrence Saul

[Y] = lle(d,K,di);

plot(Y(1,:), Y(2,:),'k.');
set(gca,'FontSize',28);
text(Y(1,:), Y(2,:),et,'FontSize',28);



 %% --- leigs function for Laplacian eigenmap.
% Written by Belkin & Niyogi, 2002.
 
[Y,V] = leigs(X, 'nn',K, di+1);

plot(Y(:,1), Y(:,2),'k.');
set(gca,'FontSize',28);

text(Y(:,1), Y(:,2),et,'FontSize',28);

%% % --- diffusionKernel function
% Written by R. Coifman & S. Lafon.
%[Y] = diffusionKernel (X,sigmaK,alpha,di)
[Y] = diffusionKernel (X,10,1,di);
%[Y] = diffusionKernel (X,2,1,di);


plot(Y(:,1), Y(:,2),'k.');
set(gca,'FontSize',28);

text(Y(:,1), Y(:,2),et,'FontSize',28);

% --- LTSA function
% Written by Zhenyue Zhang & Hongyuan Zha, 2004.
% Reference: http://epubs.siam.org/sam-bin/dbq/article/41915
 [Y,NI] = ltsa(X,di,K);

 plot(Y(1,:), Y(2,:),'k.');
set(gca,'FontSize',28);

text(Y(1,:), Y(2,:),et,'FontSize',28);
