%%%  Diffusion maps swiss roll example

% Read the Swiss Roll data set (2000 points in R^3)
filename="sr2000x3.dat"   % 2000 rows, 3 coordinates per row
separator='';       % default separator is whitespace
sr=dlmread(filename,separator);  % be in the correct folder!
[nrow,ncol]=size(sr)  % should be 2000 rows, 3 columns

% Display the data in a (blue) 3D scatterplot
plot3(sr(:,1),sr(:,2),sr(:,3),"b.");

% Take, perhaps, a subset of the data to save time
n=400;  % set n=nrow to use all the data
V=sr(1:n,:);
% Display the selected data in a (red) 3D scatterplot
plot3(V(:,1),V(:,2),V(:,3),"r.");


% Find the Euclidean distances between all distinct pairs
edists=zeros(1,n*(n-1)/2); % upper triangle
m=1;
for i=1:(n-1)
 for j=(i+1):n
  edists(m)=norm(V(i,:)-V(j,:));
  m=m+1;
 end
end
edists=sort(edists);  % increasing order
figure; plot(edists)  % for visualization


% Find the Euclidean distances to nearest neighbors
ennds=zeros(n,n); % full nxn matrix
for i=1:n
 for j=1:n
  ennds(i,j)=norm(V(i,:)-V(j,:));
 end
 ennds(i,:)=sort(ennds(i,:));
end
figure; plot(min(ennds)); hold on; plot(max(ennds))  % min,max to kth nearest neighbor


% Gaussian kernel
sig = max(ennds)(1+10)  % distance to 10th neighbor
kernel=@(v1,v2) exp(-norm(v1-v2)^2 /sig^2); % kernel function

% Kernel matrix
K=zeros(n,n);
for i=1:n
 for j=1:n
   K(i,j)= kernel(V(i,:),V(j,:));
 end
end

% Degree matrix
D=zeros(n,n);
for i=1:n
  D(i,i)=sum(K(i,:));
end

% Transition matrix
P=inv(D)*K;
figure; imagesc(P);  % visualize

% Stationary distribution
pd= diag(D)' ;  % diagonal part of matrix D, transpose to row vector
pd= pd/sum(pd); % normalize
sum(pd)        % should be 1: pd is a pdf
norm(pd*P-pd)  % should be 0: pd is stationary for P

% Symmetrize
sPd=diag(sqrt(pd)); isPd=diag(sqrt(1./pd));
A=sPd*P*isPd;
norm(A-A')  % should be zero, since A is symmetric

% Find first few eigenvectors Theta and eigenvalues Lambda of A
[Theta,Lambda]=eigs(A);

% Find eigenvectors for P
Psi = isPd*Theta;

% Diffusion map Psi_k(v_i) is the ith row
% Plot in 2D at various scales k
x=2; y=3; %  first two nonconstant columns
k=1;     Psik=(Psi*Lambda^1); figure; plot(Psik(:,x),Psik(:,y),".");
k=16;    Psik=(Psi*Lambda^k); figure; plot(Psik(:,x),Psik(:,y),".");
k=64;    Psik=(Psi*Lambda^k); figure; plot(Psik(:,x),Psik(:,y),".");
k=256;   Psik=(Psi*Lambda^k); figure; plot(Psik(:,x),Psik(:,y),".");
k=1024;  Psik=(Psi*Lambda^k); figure; plot(Psik(:,x),Psik(:,y),".");
k=4096;  Psik=(Psi*Lambda^k); figure; plot(Psik(:,x),Psik(:,y),".");

% Plot in 3D at various scales
x=2; y=3; z=4;
k=1;     Psik=(Psi*Lambda^k); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z),".");
k=16;    Psik=(Psi*Lambda^k); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z),".");
k=64;    Psik=(Psi*Lambda^k); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z),".");
k=256;   Psik=(Psi*Lambda^k); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z),".");
k=1024;  Psik=(Psi*Lambda^k); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z),".");
k=4096;  Psik=(Psi*Lambda^k); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z),".");


% Display powers of the transition matrix as images
M=P;  % initialize a new matrix to raise to dyadic powers
imagesc(M);  % M as an image, scaled to fill the color map
M=M*M; imagesc(M);  % image of M^2
M=M*M; imagesc(M);  % image of M^4
M=M*M; imagesc(M);  % image of M^8
M=M*M; imagesc(M);  % ...repeat this to see what happens

%% The following code will be slow.
% Implement the diffusion map with gaussian kernel
% as a function dmgk(V,sig,k) and plot the (partial)
% swiss roll in 2D at fixed k=32 and values sig=1,2,4,8:
x=2; y=3; %  first two nonconstant columns
Psik=dmgk(V,1,32); figure; plot(Psik(:,x),Psik(:,y),".");
Psik=dmgk(V,2,32); figure; plot(Psik(:,x),Psik(:,y),".");
Psik=dmgk(V,4,32); figure; plot(Psik(:,x),Psik(:,y),".");
Psik=dmgk(V,8,32); figure; plot(Psik(:,x),Psik(:,y),".");