% GraphTook - 2007-2012 Glotin & Chouchane
% Dyni LSIS USTV - {glotin, mathieu.chouchane} @ univ-tln.fr

% Détrmination du nombre de colorations possibles d'un graphe avec L
% couleurs (polynome chromatique).
% Par exemple, le graphe de 6 sommets ci-dessous, est coloriable avec L
% couleurs de n facons differentes, n donne par ce polynome chromatique :
% P(L) = [L]3 + 3*[L]4 + [L]5, avec [L]N le nombre de colorations a L
% couleurs de la clique KN, [L]N = L! / (L-N)! .

% Ce graphe est donc coloriable de 6 facons avec 3 couleurs, et de 96
% facons avec 4 couleurs,... Les valeurs de P pour L de 1 a 10 sont donnees
% dans la courbe ci-contre.

% Ci-dessous nous proposons une implementation par double recursion de
% l'algorithme de Hajos pour determiner ce polynome chromatique. Ce code
% (programmation vectorielle) est executable en Octave (ou MATLAB).
% Le calcul est rapide, meme pour des graphes de quelque dizaines de
% sommets. Des heuristiques d'acceleration sont discutees en fin de module.


function [pcc, pm] = polychro(G, pfm)
% le graphe est alors coloriable avec ncoul en pcc(ncoul) facons
% (ca peut etre 0 facon).
% pm est la profondeur max atteinte dans la decomposition de Hajos
% evaluer_polychro.m permet de calculer pcc(ncoul)
% pcc est construit recursivement, a chqe dim N, il donne C(N)*LAMBDA_N
% (C) Glotin dec 2009
if nargin == 1
    pfm = 0;
    % profondeur de l'arbre de Hajos pour calculer ce polynome chromatique
end
c = clique(G); % renvoi 0 si G n'est pas un clique, sinon N
if (c ~= 0)
    pcc = zeros(1,20); % on suppose que le degre max du polynome est 20
    pcc(c) = pcc(c)+1;
    pm = pfm;
else
    [i, j] = paire_disjointe(G);
    [pccg, pmg]= polychro(relie(G,i,j), pfm+1);
    [pccd, pmd]= polychro(contracte(G,i,j), pfm+1);
    pm = max(pmg, pmd);
    pcc = pccg + pccd;
end % if
end % funct

%%%%%%%%%%%%%%%
function C = contracte(G,i,j)
% matrice adj de G (simple sans boucle)
% renvoie C, contracte de G en (i,j)
% remarque : G(i,j) = 1 alors bouclage
C = G;
C(j,:) = or(C(i,:),C(j,:));
C(i,:) = [];
C(:,j) = or(C(:,i),C(:,j));
C(:,i) = [];
end % funct


function R = relie(G,i,j)
% matrice d'adjacence de G (simple sans boucle)
% renvoie R, le relie de G en (i,j)
R = G;
R(i,j) = 1;
R(j,i) = 1;
end % funct


function [x, y] = paire_disjointe(G)
% le rand est pour tirer au hasard la paire de sommets qui sont disjoints
[I, J] = find(G == 0);
n = length(I);
A = randperm(n);
k = 1;
while I(A(k)) == J(A(k))
    k = k + 1;
end % while
x = I(A(k));
y = J(A(k));
end %funct


function n = clique(G)
% matrice d'adjacence de G (simple sans boucle)
% renvoie 0 si G n'est pas une clique
% sinon renvoie n car G = Kn
n = size(G,1) ;
if sum(G(:)) ~= n * (n-1) % sum(X(:))
    n = 0;
end
end % funct


function y = lambdaN(L,N)
% Donne [Lambda]N facon(s) de colorier KN en lambda couleur(s)
% L = le nbre de couleur, N
if(L < N)
    y = 0;
elseif(L == N)
    y = prod(1:L);
else
    y = prod(1:L) / prod(1:L-N);
end
end % funct
