%SYMMETRIC_NORM computes the symmetric normalization of a nonnegative matrix

function OUT = symmetric_norm(W)

%INPUT: W = nonnegative matrix 

%OUT = symmetric-normalized W  

% NOTE: If W(i,:)is identically zero, then so is OUT(i,:)


[r,c] = size(W);

if r ~= c
   
   disp('matrix not square')
   
   break
   
end


W = (W + W')/2 ;

z = sum(W) ;

I = find(z ~= 0);

rr = length(I);

WW = full(W(I,I));  %Eliminate zero rows and columns from W

S = sqrtm(inv(diag(sum(WW))));

SS = zeros(r,c);

SS(I,I) = S;

OUT = SS*W*SS;