%%% Commands to compute diffusion maps

% Generate a data set V
dt=0.2; t=0:dt:(200*dt);  
Vx=cos(t); Vy=sin(t); Vz=t/(2*pi);
plot3(Vx,Vy,Vz);  % display the result
V=[Vx' Vy' Vz'];  % each row is a point in R3
n=length(Vx);   % number of data points
% Check distances: adjacent << after-one-turn?
norm(V(1,:)-V(2,:))  % adjacent points on spiral
norm(V(1,:)-V(1+round(2*pi/dt),:))  % after one turn

% Gaussian kernel
sig = 0.2;  % approximately the distance between neighbors
kernel=@(v1,v2) exp(-norm(v1-v2)^2 /sig^2); % kernel function
K=zeros(n,n);  % kernel matrix
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;

% 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 eigenvectors Theta and eigenvalues Lambda of A
[Theta,Lambda]=eigs(A);  % default: top 6 eigenvalues, nx6 Theta, 6x6 Lambda

% Find eigenvectors for P
Psi = isPd*Theta;

% Diffusion map Psi_k(v_i) is the ith row:
k=4; Psik=(Psi*Lambda^k);

% Plot the first two nonconstant columns:
dfx=Psik(:,2); dfy=Psik(:,3);
plot( dfx, dfy );

% Plot in 2D at various scales
x=2; y=3;
Psik=(Psi*Lambda^1); figure; plot(Psik(:,x),Psik(:,y));
Psik=(Psi*Lambda^16); figure; plot(Psik(:,x),Psik(:,y));
Psik=(Psi*Lambda^64); figure; plot(Psik(:,x),Psik(:,y));
Psik=(Psi*Lambda^256); figure; plot(Psik(:,x),Psik(:,y));
Psik=(Psi*Lambda^4096); figure; plot(Psik(:,x),Psik(:,y));
Psik=(Psi*Lambda^32768); figure; plot(Psik(:,x),Psik(:,y));

% Plot in 3D at various scales
x=2; y=3; z=4;
Psik=(Psi*Lambda^1); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z));
Psik=(Psi*Lambda^16); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z));
Psik=(Psi*Lambda^64); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z));
Psik=(Psi*Lambda^256); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z));
Psik=(Psi*Lambda^4096); figure; plot3(Psik(:,x),Psik(:,y),Psik(:,z));
Psik=(Psi*Lambda^32768); 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