chromatic_poly

 chrompoly(g) --- find the chromatic polynomial of g
 Warning: This algorithm is slow and unusable except for small graphs.
 Author: James Preen

0001 function out = chromatic_poly(g)
0002 % chrompoly(g) --- find the chromatic polynomial of g
0003 % Warning: This algorithm is slow and unusable except for small graphs.
0004 % Author: James Preen
0005 
0006 out = chrompoly(double(matrix(g)));
0007 
0008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0009 % Here is James Preen's .m file very lightly edited %
0010 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0011 function out = chrompoly(A)
0012 
0013 out=[];
0014 n=size(A,1);
0015 e=sum(sum(A))/2;
0016 
0017 if e==0,
0018     out(n+1)=0;out(1)=1;
0019     %  disp(['that is k bar ' num2str(n)]);
0020 else
0021     somenew=1;
0022     nb=ceil(rand(1)*n);
0023     while somenew,
0024         nb2=rem(find(A(nb',:)')-1,n)+1;
0025         nb2=nb2';
0026         nb3=unique([nb nb2]);
0027         if size(nb3,2)==size(nb,2), somenew=0;end;
0028         nb=nb3;
0029     end
0030     if size(nb,2)<n,
0031         %       disp(['disconnected set: ' num2str([nb])]);
0032         nb2=setdiff(1:n,nb);
0033         A1=A;A1(:,nb2)=[];A1(nb2,:)=[];
0034         A2=A;A2(:,nb)=[];A2(nb,:)=[];
0035         out=conv(chrompoly(A1),chrompoly(A2));
0036     else
0037         valseq=sum(A);
0038         % v1 = min(find(valseq==max(valseq)));
0039         v1 = find(valseq==max(valseq),1);
0040         if valseq(v1)==n-1,
0041             %   disp(['deleting vertex ' num2str([v1])]);
0042             A(v1,:)=[];A(:,v1)=[];
0043             p=chrompoly(A);
0044             bp=[1 -1];tp=1;
0045             len=size(p,2);m3a=[];
0046             for i=1:size(p,2)
0047                 m3a=[0 m3a] + p(len+1-i)*tp;
0048                 tp=conv(tp,bp);
0049             end;
0050             out=[m3a 0];
0051         else
0052     
0053             if e < n*(n-1)/2,
0054                 %     %choose 2 random vertices
0055                 %     p=find(valseq>0);
0056                 %     v1=p(ceil(rand(1)*size(p,2)));
0057                 %     nb1=find(A(v1,:)==1);
0058                 %     v2=nb1(ceil(rand(1)*valseq(v1)));
0059 
0060                 nb1=find(A(v1,:)==1);
0061                 nbv=valseq(nb1);
0062                 v2=nb1(min(find(nbv==max(nbv))));
0063 
0064                 %deletion
0065                 %   disp(['deleting edge ' num2str([v1 v2])]);
0066                 A1=A;
0067                 A1(v1,v2)=0;A1(v2,v1)=0;
0068                 m1=chrompoly(A1);
0069     
0070                 %contraction
0071                 %    disp(['contracting edge ' num2str([v1 v2])]);
0072                 A2=A;
0073                 A2(v1,:)=(A2(v2,:)+A2(v1,:))>0;
0074                 A2(:,v1)=(A2(:,v2)+A2(:,v1))>0;
0075                 A2(v1,v1)=0;A2(:,v2)=[];A2(v2,:)=[];
0076                 m2=chrompoly(A2);
0077     
0078                 s1=size(m1,2);
0079                 s2=size(m2,2);
0080                 pd=zeros(1,abs(s2-s1));
0081     
0082                 if s1>s2,
0083                     out=m1 - [pd m2];
0084                 elseif s2>s1,
0085                     out=[pd m1] - m2;
0086                 else
0087                     out= m1 - m2;
0088                 end
0089             else
0090     
0091                 nb1=find(A(v1,:)==0);
0092                 nbv=valseq(nb1);
0093                 v2=v1;
0094                 while v2==v1, v2=nb1(min(find(nbv==max(nbv)))); end;
0095 
0096                 %deletion
0097                 %   disp(['adding edge ' num2str([v1 v2])]);
0098                 A1=A;
0099                 A1(v1,v2)=1;A1(v2,v1)=1;
0100                 m1=chrompoly(A1);
0101     
0102                 %contraction
0103                 %    disp(['contracting edge ' num2str([v1 v2])]);
0104                 A2=A;
0105                 A2(v1,:)=(A2(v2,:)+A2(v1,:))>0;
0106                 A2(:,v1)=(A2(:,v2)+A2(:,v1))>0;
0107                 A2(v1,v1)=0;A2(:,v2)=[];A2(v2,:)=[];
0108                 m2=chrompoly(A2);
0109     
0110                 s1=size(m1,2);
0111                 s2=size(m2,2);
0112                 pd=zeros(1,abs(s2-s1));
0113     
0114                 if s1>s2,
0115                     out=m1 + [pd m2];
0116                 elseif s2>s1,
0117                     out=[pd m1] + m2;
0118                 else
0119                     out= m1 + m2;
0120                 end
0121 
0122             end; %del-cont  / add-id
0123         end;
0124     end;
0125 end; %if e==0
0126 
0127 
0128 
0129 %%% OLD VERSION
0130 %
0131 % n = nv(g);
0132 % m = ne(g);
0133 %
0134 % if m==0
0135 %     p = zeros(1,n+1);
0136 %     p(1) = 1;
0137 %     return
0138 % end
0139 %
0140 % elist = edges(g);
0141 % u = elist(1,1);
0142 % v = elist(1,2);
0143 %
0144 % h = graph;
0145 % copy(h,g);
0146 %
0147 % delete(h,u,v);
0148 % p1 = chrompoly(h);
0149 % contract(h,u,v);
0150 % p2 = chrompoly(h);
0151 %
0152 % p = p1 - [0,p2];
0153 %
0154 % free(h);