iso

 [yn,p] = iso(g,h,options) --- is g isomorphic to h?
 Given graphs g and h, determine whether or not the graphs are isomorphic, 
 and if so, return a permutation p such that renumber(g,p) makes g==h.
 
 Returns yn = 1 if isomorphic and yn = 0 if not. The optional p is a
 permutation giving the isomorphism.

 The optional third argument, options, is used to set various parameters
 guiding the search for an isomorphism. Here are the fields that can be
 specified in options and their default values.

 options.verbose --- turn off (0, default) or on (1) verbose output.

 options.eig --- cospectrality tests. Set this to...
     negative value: skip this test
     zero: find characteristic polynomial 
     positive: calculate eigenvalues and use this value for a threshold
     for equality.   default value = nv * 20 * eps

 options.pretest --- do basic tests such as number of vertices, edges.
     Set to 1 for on (default) and 0 for off (but don't do that!); these
     are generally very fast
 
 options.components --- check that the number and size of components are
     the same. 1 = check, 0 = don't check. Default is 0

 options.distance --- use the distance matrix to distinguish vertex types.
     1 = check (default), 0 = don't check. Default is 1. This can speed up
     isomorphism testing for regular graphs.

0001 function [yn,p] = iso(g,h,options)
0002 % [yn,p] = iso(g,h,options) --- is g isomorphic to h?
0003 % Given graphs g and h, determine whether or not the graphs are isomorphic,
0004 % and if so, return a permutation p such that renumber(g,p) makes g==h.
0005 %
0006 % Returns yn = 1 if isomorphic and yn = 0 if not. The optional p is a
0007 % permutation giving the isomorphism.
0008 %
0009 % The optional third argument, options, is used to set various parameters
0010 % guiding the search for an isomorphism. Here are the fields that can be
0011 % specified in options and their default values.
0012 %
0013 % options.verbose --- turn off (0, default) or on (1) verbose output.
0014 %
0015 % options.eig --- cospectrality tests. Set this to...
0016 %     negative value: skip this test
0017 %     zero: find characteristic polynomial
0018 %     positive: calculate eigenvalues and use this value for a threshold
0019 %     for equality.   default value = nv * 20 * eps
0020 %
0021 % options.pretest --- do basic tests such as number of vertices, edges.
0022 %     Set to 1 for on (default) and 0 for off (but don't do that!); these
0023 %     are generally very fast
0024 %
0025 % options.components --- check that the number and size of components are
0026 %     the same. 1 = check, 0 = don't check. Default is 0
0027 %
0028 % options.distance --- use the distance matrix to distinguish vertex types.
0029 %     1 = check (default), 0 = don't check. Default is 1. This can speed up
0030 %     isomorphism testing for regular graphs.
0031 %
0032 
0033 %% Check the options structure and set default values
0034 
0035 if nargin < 3
0036     options.verbose = 0;
0037 end
0038 
0039 if ~isfield(options,'verbose')
0040     options.verbose = 0;
0041 end
0042 
0043 if ~isfield(options,'eig')
0044     options.eig = nv(g) * 20 *eps;
0045 end
0046 
0047 if ~isfield(options,'pretest')
0048     options.pretest = 1;
0049 end
0050 
0051 if ~isfield(options,'components')
0052     options.components = 0;
0053 end
0054 
0055 if ~isfield(options,'distance')
0056     options.distance = 1;
0057 end
0058 
0059 if ~isfield(options,'debug')
0060     options.debug = false;
0061 end
0062 
0063 
0064 verb = options.verbose;
0065 
0066 if verb
0067     disp('Searching for an isomorphism')
0068 end
0069 
0070 % default in case the graphs are not isomorphic
0071 yn = false;
0072 p = permutation(0);
0073 
0074 %%
0075 % We begin with a set of basic checks that can easily determine if the
0076 % graphs are not isomorphic.
0077 
0078 
0079 % check if the graphs are simply equal
0080 
0081 if verb
0082     disp('Checking if the graphs are identical')
0083     tic
0084 end
0085 
0086 if (g==h)
0087     yn = true;
0088     p = permutation(nv(g));
0089     if verb
0090         toc
0091     end
0092     return
0093 end
0094 
0095 
0096 % number of vertices & edges
0097 if options.pretest
0098     if ~iso_pretest(g,h,options)
0099         if verb
0100             toc
0101         end
0102         return
0103     end
0104 end
0105 
0106 
0107 %% Create vertex distinguishing information. This is kept in n-row matrices
0108 % Mg and Mh.
0109 
0110 % first set of distinguishing mark is the degree sequence of its neighbors
0111 
0112 if verb
0113     disp('Calculating vertex neighbor degree sequences')
0114 end
0115 
0116 n = nv(g);
0117 
0118 Mg = deg(g);
0119 Mg = Mg(:);
0120 Mh = deg(h);
0121 Mh = Mh(:);
0122     
0123 Mg = refine_types(g,Mg);
0124 Mh = refine_types(h,Mh);
0125 
0126 % check if they're the same
0127 
0128 if verb
0129     toc
0130 end
0131 
0132 if any(any( sortrows(Mg) ~= sortrows(Mh) ))
0133     return
0134 end
0135 
0136 % the second set of distinguishing marks is from the distance matrix
0137 
0138 if options.distance
0139     if verb
0140         disp('Calculating distance matrices')
0141     end
0142     Dg = sort(dist(g),2);
0143     Dh = sort(dist(h),2);
0144 
0145     % append to Mg and Mh
0146 
0147     Mg = [Mg, Dg(:,2:end)];
0148     Mh = [Mh, Dh(:,2:end)];
0149     
0150     if verb
0151         toc
0152     end
0153     
0154 end
0155 
0156 
0157 % make sure they're the same
0158 if any(any( sortrows(Mg) ~= sortrows(Mh) ))
0159     if options.verbose
0160         disp('Graphs have different distance structure')
0161     end
0162     return
0163 end
0164 
0165 
0166 %% Now we assign a "type" to each vertex of the graph based on its row in
0167 % Mg (or Mh). Vertices of different types can NOT be linked by an
0168 % isomorphism.
0169 
0170 
0171 % determine type number for each vertex
0172 
0173 if verb
0174     disp('Refining')
0175 end
0176 
0177 gtypes = matrix2types(Mg);
0178 htypes = matrix2types(Mh);
0179 
0180 gtypes = ultra_refine_types(g,gtypes);
0181 htypes = ultra_refine_types(h,htypes);
0182 
0183 
0184 ntypes = max(gtypes);
0185 
0186 if any(sort(gtypes) ~= sort(htypes))
0187     yn = 0;
0188     if verb 
0189         toc
0190     end
0191     return
0192 end
0193 
0194 
0195 if verb
0196     disp(['Number of different types of vertices = ', int2str(ntypes)])
0197     toc
0198 end
0199 
0200 [new_gtypes,new_htypes] = iso_match(g,h,gtypes,htypes,options);
0201 
0202 if max(new_gtypes) < n || max(new_htypes) < n
0203     return
0204 end
0205 
0206 p1 = permutation(new_gtypes);
0207 p2 = permutation(new_htypes);
0208 
0209 p = inv(p2)*p1;
0210 
0211 gcopy = graph;
0212 copy(gcopy, g);
0213 renumber(gcopy,p);
0214 if gcopy==h
0215     yn = 1;
0216     if options.verbose
0217         disp('Graphs are isomorphic');
0218         toc
0219     end
0220 else
0221     yn = 0;
0222     p = permutation(0);
0223     if options.verbose
0224         disp('Graphs are NOT isomorphic');
0225         toc
0226     end
0227 end
0228 free(gcopy);
0229 
0230 
0231 %% main work engine
0232 function [new_gtypes, new_htypes] = iso_match(g,h,gtypes,htypes,options)
0233 
0234 ntypes = max(gtypes); 
0235 n = nv(g);
0236 
0237 if ntypes == n
0238     new_gtypes = gtypes;
0239     new_htypes = htypes;
0240     return
0241 end
0242 
0243 % find a g-vertex of the most frequent type
0244 
0245 type_count = zeros(1,ntypes);
0246 
0247 % find a vertex of the most numerous type
0248 for k=1:ntypes
0249     type_count(k) = length(find(k==gtypes));
0250 end
0251 max_idx = find(type_count == max(type_count));
0252 max_type = max_idx(1);
0253 
0254 % get a g-vertex of type max_type; call it v
0255     
0256 v = find(gtypes == max_type);
0257 v = v(1);
0258 
0259 new_num = ntypes + 1; % create a new type for v
0260 gtypes(v) = new_num;
0261 new_gtypes = ultra_refine_types(g,gtypes);
0262 
0263 
0264 % find all h vertices of the same type as v
0265 % call that list wlist
0266 
0267 wlist = find(htypes == max_type);
0268 
0269 nw = length(wlist);
0270 for i = 1:nw
0271     % check if G-v is plausible match to H-w
0272     w = wlist(i);
0273     
0274     if options.verbose
0275         disp(['Trying ', int2str(v), ' --> ', int2str(w)])
0276         toc
0277     end
0278     
0279     if options.pretest
0280         g_v = graph; copy(g_v,g); delete(g_v,v);
0281         h_w = graph; copy(h_w,h); delete(h_w,w);
0282         if (~iso_pretest(g_v,h_w,options))
0283             if options.verbose
0284                 disp('failed pretest')
0285             end
0286             free(g_v);
0287             free(h_w);
0288             continue        
0289         end
0290         free(g_v);
0291         free(h_w);
0292     
0293         if options.verbose
0294             disp('Pretest passed')
0295         end
0296     end
0297 
0298     % plausible match, so assign new type to w and recurse
0299     htypes_tmp = htypes;
0300     htypes_tmp(w) = new_num;
0301     htypes_tmp = ultra_refine_types(h,htypes_tmp);
0302     
0303     
0304     % htypes_tmp should be same as gtypes (up to sort order)
0305     if any(sort(htypes_tmp) ~= sort(new_gtypes))
0306         if options.verbose
0307             disp('Attempted match failed')
0308         end
0309         continue;
0310     end
0311     
0312     
0313     if options.verbose
0314         disp(['Number of vertex types is now ', int2str(max(htypes_tmp))])
0315         toc
0316     end
0317     
0318     % recurse
0319     [new_gtypes, new_htypes] = iso_match(g,h,new_gtypes,htypes_tmp,options);
0320     
0321     if (max(new_gtypes) == n) && (max(new_htypes)==n)
0322         p1 = permutation(new_gtypes);
0323         p2 = permutation(new_htypes);
0324 
0325         p = inv(p2)*p1;
0326 
0327         gcopy = graph;
0328         copy(gcopy, g);
0329         renumber(gcopy,p);
0330         if gcopy==h
0331             yn = 1;
0332         else
0333             yn = 0;
0334         end
0335         free(gcopy);
0336         
0337         if yn
0338             return
0339         end
0340     end
0341 end
0342 
0343 new_gtypes = gtypes;
0344 new_htypes = htypes;
0345     
0346 
0347 %% matrix2types -- convert rows of M to type numbers
0348 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0349 
0350 function tlist = matrix2types(M)
0351 % convert the rows of M into a column of row types
0352 
0353 T = unique(sortrows(M),'rows');
0354 
0355 n = size(M,1);
0356 
0357 % determine type number for each row
0358 
0359 tlist = zeros(n,1);
0360 
0361 for v=1:n
0362     r = M(v,:);
0363     tlist(v) = find_row(r,T);
0364 end
0365 
0366 
0367 
0368 %% refine_types
0369 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0370 
0371 function newtypes = refine_types(g,types)
0372 % newtypes = refine_types(g,types) -- refine types based on neighbors' types
0373 % given an n-vector of vertex types, augment this by considering the types
0374 % of the neighbors
0375 
0376 
0377 maxd = max(deg(g));
0378 n = nv(g);
0379 
0380 M = zeros(n,maxd+1);
0381 
0382 for v=1:n
0383     Nv = neighbors(g,v);
0384     row = types(Nv);
0385     row = sort(row);
0386     row = row(:)';  % make sure it's a row vector
0387     row = [types(v), row, zeros(1, maxd-length(row))];
0388     M(v,:) = row;
0389 end
0390 
0391 
0392 newtypes = matrix2types(M);
0393 
0394 %% ultra_refine_types -- repeatedly apply refine_types
0395 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0396 function newtypes = ultra_refine_types(g,types)
0397 % repeatedly apply refine_types until no change
0398 
0399 while true
0400     newtypes = refine_types(g,types);
0401     if all(newtypes == types)
0402         return
0403     end
0404     types = newtypes;
0405 end
0406 
0407 
0408 
0409 %% iso_pretest
0410 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0411 function yn = iso_pretest(g,h,options)
0412 % basic tests that can rule out isomorphism
0413 
0414 % # of vertices and edges
0415 
0416 yn = 0;
0417 if any(size(g) ~= size(h))
0418     return
0419 end
0420 
0421 % degree sequences
0422 
0423 dg = sort(deg(g));
0424 dh = sort(deg(h));
0425 if any(dg ~= dh)
0426     if options.verbose
0427         disp('Graphs have different degree sequences')
0428     end
0429     return
0430 end  
0431 
0432 % components
0433 
0434 if options.components
0435 
0436 
0437     gc = components(g);
0438     hc = components(h);
0439 
0440     if np(gc) ~= np(hc)
0441         return
0442     end
0443 
0444     %  components must all have the same sizes
0445     gcp = parts(gc);
0446     hcp = parts(hc);
0447 
0448     gx = zeros(1,np(gc));
0449     gy = gx;
0450 
0451     for k=1:np(gc)
0452         gx(k) = length(gcp{k});
0453         gy(k) = length(hcp{k});
0454     end
0455     gx = sort(gx);
0456     gy = sort(gy);
0457     if any(gx ~= gy)
0458         if options.verbose
0459             disp('Graphs have different component sizes')
0460         end
0461         return
0462     end
0463 end % end component subtest
0464 
0465 % cospectrality tests
0466 
0467 if options.eig == 0
0468     
0469     A = double(matrix(g));
0470     B = double(matrix(h));
0471     pA = round(poly(A));
0472     pB = round(poly(B));
0473     if any(pA ~= pB)
0474         if options.verbose
0475             disp('Graphs have different characteristic polynomials')
0476         end
0477         return;
0478     end
0479 end
0480 
0481 
0482 if options.eig > 0
0483     eA = sort(eig(g));
0484     eB = sort(eig(h));
0485     diff = norm(eA-eB);
0486     if (diff > options.eig)
0487         if options.verbose
0488             disp(['Graphs have different spectra (tol = ', ... 
0489                 num2str(options.eig),')']);
0490         end
0491         return
0492     end
0493 end
0494 
0495 
0496 yn = 1;