0001 function [yn,p] = iso(g,h,options)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
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
0071 yn = false;
0072 p = permutation(0);
0073
0074
0075
0076
0077
0078
0079
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
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
0108
0109
0110
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
0127
0128 if verb
0129 toc
0130 end
0131
0132 if any(any( sortrows(Mg) ~= sortrows(Mh) ))
0133 return
0134 end
0135
0136
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
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
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
0167
0168
0169
0170
0171
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
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
0244
0245 type_count = zeros(1,ntypes);
0246
0247
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
0255
0256 v = find(gtypes == max_type);
0257 v = v(1);
0258
0259 new_num = ntypes + 1;
0260 gtypes(v) = new_num;
0261 new_gtypes = ultra_refine_types(g,gtypes);
0262
0263
0264
0265
0266
0267 wlist = find(htypes == max_type);
0268
0269 nw = length(wlist);
0270 for i = 1:nw
0271
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
0297 htypes_tmp = htypes;
0298 htypes_tmp(w) = new_num;
0299 htypes_tmp = ultra_refine_types(h,htypes_tmp);
0300
0301
0302
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
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
0346
0347
0348 function tlist = matrix2types(M)
0349
0350
0351 T = unique(sortrows(M),'rows');
0352
0353 n = size(M,1);
0354
0355
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
0367
0368
0369 function newtypes = refine_types(g,types)
0370
0371
0372
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(:)';
0385 row = [types(v), row, zeros(1, maxd-length(row))];
0386 M(v,:) = row;
0387 end
0388
0389
0390 newtypes = matrix2types(M);
0391
0392
0393
0394 function newtypes = ultra_refine_types(g,types)
0395
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
0408
0409 function yn = iso_pretest(g,h,options)
0410
0411
0412
0413
0414 yn = 0;
0415 if any(size(g) ~= size(h))
0416 return
0417 end
0418
0419
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
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
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
0462
0463
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;