Home > matgraph > @graph > iso.m

iso

PURPOSE ^

[yn,p] = iso(g,h,options) --- is g isomorphic to h?

SYNOPSIS ^

function [yn,p] = iso(g,h,options)

DESCRIPTION ^

 [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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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             continue
0287         end
0288         free(g_v);
0289         free(h_w);
0290     
0291         if options.verbose
0292             disp('Pretest passed')
0293         end
0294     end
0295 
0296     % plausible match, so assign new type to w and recurse
0297     htypes_tmp = htypes;
0298     htypes_tmp(w) = new_num;
0299     htypes_tmp = ultra_refine_types(h,htypes_tmp);
0300     
0301     
0302     % htypes_tmp should be same as gtypes (up to sort order)
0303     if any(sort(htypes_tmp) ~= sort(new_gtypes))
0304         if options.verbose
0305             disp('Attempted match failed')
0306         end
0307         continue;
0308     end
0309     
0310     
0311     if options.verbose
0312         disp(['Number of vertex types is now ', int2str(max(htypes_tmp))])
0313         toc
0314     end
0315     
0316     % recurse
0317     [new_gtypes, new_htypes] = iso_match(g,h,new_gtypes,htypes_tmp,options);
0318     
0319     if (max(new_gtypes) == n) && (max(new_htypes)==n)
0320         p1 = permutation(new_gtypes);
0321         p2 = permutation(new_htypes);
0322 
0323         p = inv(p2)*p1;
0324 
0325         gcopy = graph;
0326         copy(gcopy, g);
0327         renumber(gcopy,p);
0328         if gcopy==h
0329             yn = 1;
0330         else
0331             yn = 0;
0332         end
0333         free(gcopy);
0334         
0335         if yn
0336             return
0337         end
0338     end
0339 end
0340 
0341 new_gtypes = gtypes;
0342 new_htypes = htypes;
0343     
0344 
0345 %% matrix2types -- convert rows of M to type numbers
0346 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0347 
0348 function tlist = matrix2types(M)
0349 % convert the rows of M into a column of row types
0350 
0351 T = unique(sortrows(M),'rows');
0352 
0353 n = size(M,1);
0354 
0355 % determine type number for each row
0356 
0357 tlist = zeros(n,1);
0358 
0359 for v=1:n
0360     r = M(v,:);
0361     tlist(v) = find_row(r,T);
0362 end
0363 
0364 
0365 
0366 %% refine_types
0367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0368 
0369 function newtypes = refine_types(g,types)
0370 % newtypes = refine_types(g,types) -- refine types based on neighbors' types
0371 % given an n-vector of vertex types, augment this by considering the types
0372 % of the neighbors
0373 
0374 
0375 maxd = max(deg(g));
0376 n = nv(g);
0377 
0378 M = zeros(n,maxd+1);
0379 
0380 for v=1:n
0381     Nv = neighbors(g,v);
0382     row = types(Nv);
0383     row = sort(row);
0384     row = row(:)';  % make sure it's a row vector
0385     row = [types(v), row, zeros(1, maxd-length(row))];
0386     M(v,:) = row;
0387 end
0388 
0389 
0390 newtypes = matrix2types(M);
0391 
0392 %% ultra_refine_types -- repeatedly apply refine_types
0393 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0394 function newtypes = ultra_refine_types(g,types)
0395 % repeatedly apply refine_types until no change
0396 
0397 while true
0398     newtypes = refine_types(g,types);
0399     if all(newtypes == types)
0400         return
0401     end
0402     types = newtypes;
0403 end
0404 
0405 
0406 
0407 %% iso_pretest
0408 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0409 function yn = iso_pretest(g,h,options)
0410 % basic tests that can rule out isomorphism
0411 
0412 % # of vertices and edges
0413 
0414 yn = 0;
0415 if any(size(g) ~= size(h))
0416     return
0417 end
0418 
0419 % degree sequences
0420 
0421 dg = sort(deg(g));
0422 dh = sort(deg(h));
0423 if any(dg ~= dh)
0424     if options.verbose
0425         disp('Graphs have different degree sequences')
0426     end
0427     return
0428 end  
0429 
0430 % components
0431 
0432 if options.components
0433 
0434 
0435     gc = components(g);
0436     hc = components(h);
0437 
0438     if np(gc) ~= np(hc)
0439         return
0440     end
0441 
0442     %  components must all have the same sizes
0443     gcp = parts(gc);
0444     hcp = parts(hc);
0445 
0446     gx = zeros(1,np(gc));
0447     gy = gx;
0448 
0449     for k=1:np(gc)
0450         gx(k) = length(gcp{k});
0451         gy(k) = length(hcp{k});
0452     end
0453     gx = sort(gx);
0454     gy = sort(gy);
0455     if any(gx ~= gy)
0456         if options.verbose
0457             disp('Graphs have different component sizes')
0458         end
0459         return
0460     end
0461 end % end component subtest
0462 
0463 % cospectrality tests
0464 
0465 if options.eig == 0
0466     
0467     A = double(matrix(g));
0468     B = double(matrix(h));
0469     pA = round(poly(A));
0470     pB = round(poly(B));
0471     if any(pA ~= pB)
0472         if options.verbose
0473             disp('Graphs have different characteristic polynomials')
0474         end
0475         return;
0476     end
0477 end
0478 
0479 
0480 if options.eig > 0
0481     
0482     A = double(matrix(g));
0483     B = double(matrix(h));
0484     eA = sort(eig(A));
0485     eB = sort(eig(B));
0486     diff = norm(eA-eB);
0487     if (diff > options.eig)
0488         if options.verbose
0489             disp(['Graphs have different spectra (tol = ', ... 
0490                 num2str(options.eig),')']);
0491         end
0492         return
0493     end
0494 end
0495 
0496 
0497 yn = 1;

Generated on Thu 13-Mar-2008 14:23:52 by m2html © 2003