function C = coloring_Dsatur(G)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% General description:
%
% This function performs heuristic graph coloring using saturation degree.
%
% See details in:
% D. Brelaz, "New methods to color vertices of a graph," Communications of
% the ACM, 22-4, pp251--256, 1979
%
% input:
%       G(n,n) - The whole graph. If G(i,j)=1, then (i,j) is an edge in G
%       preordered - A subgraph that has already been ordered by GLF.
%
% output: 
%       C(n,1) - Color array. C(i) is the color of i-th vertex in the
%     
% Written by: Wei Wang
% Machine Learning for Signal Processing Lab
% University of Maryland, Baltimore County
% 8/21/2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=size(G,1);
if n==0 % empty graph
    C=[];
    return;
end

if n==1 % single vertex graph
    C=1;
    return;
end

C=zeros(n,1); % colors of all vertcies
deg=sum(G);
[v_deg_des,v_deg_des_id]=sort(deg,'descend'); % sort vertices with degree descending

K=max(sum(G))+1; % A safest estimate of the largest number of saturated colors 
Gsat=zeros(n,K); % Saturated colors
Dsat=zeros(n,1); % Degree of saturation

% Color the largest degree vertex with color 1
v= v_deg_des_id(1);
C(v)=1;
nbs=find(G(v,:)>0); %neighbour of 'v'
Dsat(nbs)=1; 
Gsat(nbs,1)=C(v); % 

% Sequentially color vertex with a maximal saturation degree by the least
% available color
while 1
    uncol=find(C==0);
    if length(uncol)==0
        break;
    end
    
    Dsatuncol=Dsat(uncol); % Saturation degree of uncolored vertices
    pos=find(Dsatuncol==max(Dsatuncol));
    v=uncol(pos(1)); % vertex with the largest saturation degree
    
    % Color v with the least available color
    vcol=Gsat(v,1:Dsat(v));
    tid=[1:K];
    tid(vcol)=[];
    col_least_available=tid(1);    
    C(v)=col_least_available;
    
    % update degree of saturation for neighbours of 'v'
    nbs=find(G(v,:)>0); %neighbour
    if numel(nbs)==0
        continue;
    end
    for i=1:length(nbs);
%        if length(find(Gsat(nbs(i),1:Dsat(nbs(i)))==C(v)))==0
        if ismember(C(v), Gsat(nbs(i),1:Dsat(nbs(i))))==0
            Dsat(nbs(i))=Dsat(nbs(i))+1;     
            tt=[Gsat(nbs(i),1:Dsat(nbs(i))-1) C(v)]; % 
            Gsat(nbs(i),1:Dsat(nbs(i)))=sort(tt,'ascend');
        end
    end
end

