function Psik=dmgk(V,sig,k)
% Compute the diffusion mapping for vertex set V
% - use the gaussian kernel function with s.d. sig 
% - use scale k
% - return the first 6 coordinates, the default from eigs()
  [n,d]=size(V);
  kernel=@(v1,v2) exp(-norm(v1-v2)^2 /sig^2); % Gaussian 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
  D=zeros(n,n);				% Degree matrix
  for i=1:n
    D(i,i)=sum(K(i,:));
  end
  P=inv(D)*K;% Transition matrix
  % Stationary distribution:
  pd= diag(D)' ;  % diagonal part of matrix D, transpose to row vector
  pd= pd/sum(pd); % normalize
  sPd=diag(sqrt(pd)); isPd=diag(sqrt(1./pd));
  A=sPd*P*isPd;   % Symmetrize
  % Find first few eigenvectors Theta and eigenvalues Lambda of A
  [Theta,Lambda]=eigs(A);
  Psi = isPd*Theta;   % eigenvectors for P
  Psik=Psi*Lambda^k;  % Diffusion map Psi_k(v_i) is the ith row
  return
end
