% example of Lsym * D^-0.5 unmixed

n=50;
eps = 3.5;
[x,y] = randomdiskpoint(n);
X = [];
x = x';
y = y';
X(1,:) = x;
X(2,:) = y;
W = ndimweightmatrix(X,eps);
D = diag(sum(W,2));

sized = size(D,1);
L = D - W ;                     % unnormalized
Lsym = eye(sized(1)) - D^-0.5 * W * D^-0.5 ;        % symmetric
Lrw = eye(sized(1)) - D^-1 * W ;

[V,lambda] = eig(Lsym); %eigenvectors for lambda = 0
V = D^-0.5 * V;

% Unmix eigenvectors
c = 0;
% first, generate E (eigenvectors with lambda=0)
E = zeros(sized,1);
for i=1:sized
    if abs(lambda(i,i)) < 0.00001
        c = c + 1;
        E(1:sized,c) = V(1:sized, i);
    end
end
% then, PA=LU factorize E and calculate imat = E*(first c rows of P*E)^-1
[L, U, P] = lu(E);
M = P*E;
M = M(1:c,1:c);
imat = E*M^-1;
imat = round(imat);

colors = [1 0 0;        % red
          0 1 0;        % light green
          0 0 1;        % dark blue
          0.95 0.95 0;  % yellow
          1 0 1;        % magenta
          0 0.95 0.95;  % light blue
          1 0.5 0;      % orange
          0.5 0 0.9;    % purple
          0.7 0.1 0.1;  % dark red
          0.1 0.5 0.1;  % dark green
          0.7 0.5 0.4;  % light brown
          0 0 0; ];     % black

% Turn indicator matrix into indicator array & color array
indic = zeros(1,sized);
col = zeros(sized,3);
for i=1:sized
    for j=1:sized
        if imat(i,j) ~= 0
            indic(i) = j;
            col(i,1:3) = colors(1+mod(j-1,length(colors)),1:3);
            break;
        end
    end
end
c;
indic;
col;

scatter(x, y, 36, col, 'filled')