1 | | """ |
2 | | N.I.C.E. - Nice (as in open source) Isomorphism Check Engine |
3 | | |
4 | | Automorphism group computation and isomorphism checking for |
5 | | graphs. |
6 | | |
7 | | This is an open source implementation of Brendan McKay's algorithm |
8 | | for graph automorphism and isomorphism. McKay released a C version |
9 | | of his algorithm, named nauty (No AUTomorphisms, Yes?) under a |
10 | | license that is not GPL compatible. Although the program is open |
11 | | source, reading the source disallows anyone from recreating |
12 | | anything similar and releasing it under the GPL. Also, many people |
13 | | have complained that the code is difficult to understand. The first |
14 | | main goal of NICE was to produce a genuinely open graph isomorphism |
15 | | program, which has been accomplished. The second goal is for this |
16 | | code to be understandable, so that computed results can be trusted |
17 | | and further derived work will be possible. |
18 | | |
19 | | To determine the isomorphism type of a graph, it is convenient to |
20 | | define a canonical label for each isomorphism class- essentially an |
21 | | equivalence class representative. Loosely (albeit incorrectly), the |
22 | | canonical label is defined by enumerating all labeled graphs, then |
23 | | picking the maximal one in each isomorphism class. The NICE |
24 | | algorithm is essentially a backtrack search. It searches through |
25 | | the rooted tree of partition nests (where each partition is |
26 | | equitable) for implicit and explicit automorphisms, and uses this |
27 | | information to eliminate large parts of the tree from further |
28 | | searching. Since the leaves of the search tree are all discrete |
29 | | ordered partitions, each one of these corresponds to an ordering of |
30 | | the vertices of the graph, i.e. another member of the isomorphism |
31 | | class. Once the algorithm has finished searching the tree, it will |
32 | | know which leaf corresponds to the canonical label. In the process, |
33 | | generators for the automorphism group are also produced. |
34 | | |
35 | | AUTHORS: |
36 | | |
37 | | - Robert L. Miller (2007-03-20): initial version |
38 | | |
39 | | - Tom Boothby (2007-03-20): help with indicator function |
40 | | |
41 | | - Robert L. Miller (2007-04-07-30): optimizations |
42 | | |
43 | | - Robert L. Miller (2007-07-07-14): PartitionStack and OrbitPartition |
44 | | |
45 | | - Tom Boothby (2007-07-14) datastructure advice |
46 | | |
47 | | - Robert L. Miller (2007-07-16-20): bug fixes |
48 | | |
49 | | REFERENCE: |
50 | | |
51 | | - [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus |
52 | | Numerantium, Vol. 30 (1981), pp. 45-87. |
53 | | |
54 | | .. note:: |
55 | | |
56 | | #. Often we assume that G is a graph on vertices |
57 | | 0,1,...,n-1. |
58 | | |
59 | | #. There is no s == loads(dumps(s)) type test since none of the |
60 | | classes defined here are meant to be instantiated for longer |
61 | | than the algorithm runs (i.e. pickling is not relevant here). |
62 | | """ |
63 | | |
64 | | #***************************************************************************** |
65 | | # Copyright (C) 2006 - 2007 Robert L. Miller <rlmillster@gmail.com> |
66 | | # |
67 | | # Distributed under the terms of the GNU General Public License (GPL) |
68 | | # http://www.gnu.org/licenses/ |
69 | | #***************************************************************************** |
70 | | |
71 | | from sage.graphs.graph import GenericGraph, Graph, DiGraph |
72 | | from sage.misc.misc import cputime |
73 | | from sage.rings.integer cimport Integer |
74 | | |
75 | | cdef class OrbitPartition: |
76 | | """ |
77 | | An OrbitPartition is simply a partition which keeps track of the |
78 | | orbits of the part of the automorphism group so far discovered. |
79 | | Essentially a union-find datastructure. |
80 | | |
81 | | EXAMPLES:: |
82 | | |
83 | | sage: from sage.graphs.graph_isom import OrbitPartition |
84 | | sage: K = OrbitPartition(20) |
85 | | sage: K.find(7) |
86 | | 7 |
87 | | sage: K.union_find(7, 12) |
88 | | sage: K.find(12) |
89 | | 7 |
90 | | sage: J = OrbitPartition(20) |
91 | | sage: J.is_finer_than(K, 20) |
92 | | True |
93 | | sage: K.is_finer_than(J, 20) |
94 | | False |
95 | | |
96 | | :: |
97 | | |
98 | | sage: from sage.graphs.graph_isom import OrbitPartition |
99 | | sage: Theta1 = OrbitPartition(10) |
100 | | sage: Theta2 = OrbitPartition(10) |
101 | | sage: Theta1.union_find(0,1) |
102 | | sage: Theta1.union_find(2,3) |
103 | | sage: Theta1.union_find(3,4) |
104 | | sage: Theta1.union_find(5,6) |
105 | | sage: Theta1.union_find(8,9) |
106 | | sage: Theta2.vee_with(Theta1, 10) |
107 | | sage: for i in range(10): |
108 | | ... print i, Theta2.find(i) |
109 | | 0 0 |
110 | | 1 0 |
111 | | 2 2 |
112 | | 3 2 |
113 | | 4 2 |
114 | | 5 5 |
115 | | 6 5 |
116 | | 7 7 |
117 | | 8 8 |
118 | | 9 8 |
119 | | """ |
120 | | |
121 | | def __new__(self, int n): |
122 | | cdef int k |
123 | | self.elements = <int *> sage_malloc( n * sizeof(int) ) |
124 | | if not self.elements: |
125 | | raise MemoryError("Error allocating memory.") |
126 | | self.sizes = <int *> sage_malloc( n * sizeof(int) ) |
127 | | if not self.sizes: |
128 | | sage_free(self.elements) |
129 | | raise MemoryError("Error allocating memory.") |
130 | | for k from 0 <= k < n: |
131 | | self.elements[k] = -1 |
132 | | self.sizes[k] = 1 |
133 | | |
134 | | def __dealloc__(self): |
135 | | sage_free(self.elements) |
136 | | sage_free(self.sizes) |
137 | | |
138 | | def find(self, x): |
139 | | """ |
140 | | Returns an element of the cell which depends only on the cell. |
141 | | |
142 | | EXAMPLE:: |
143 | | |
144 | | sage: from sage.graphs.graph_isom import OrbitPartition |
145 | | sage: K = OrbitPartition(20) |
146 | | |
147 | | 0 and 1 begin in different cells:: |
148 | | |
149 | | sage: K.find(0) |
150 | | 0 |
151 | | sage: K.find(1) |
152 | | 1 |
153 | | |
154 | | Now we put them in the same cell:: |
155 | | |
156 | | sage: K.union_find(0,1) |
157 | | sage: K.find(0) |
158 | | 0 |
159 | | sage: K.find(1) |
160 | | 0 |
161 | | """ |
162 | | return self._find(x) |
163 | | |
164 | | cdef int _find(self, int x): |
165 | | if self.elements[x] == -1: |
166 | | return x |
167 | | self.elements[x] = self._find(self.elements[x]) |
168 | | return self.elements[x] |
169 | | |
170 | | def union_find(self, a, b): |
171 | | """ |
172 | | Merges the cells containing a and b. |
173 | | |
174 | | EXAMPLE:: |
175 | | |
176 | | sage: from sage.graphs.graph_isom import OrbitPartition |
177 | | sage: K = OrbitPartition(20) |
178 | | |
179 | | 0 and 1 begin in different cells:: |
180 | | |
181 | | sage: K.find(0) |
182 | | 0 |
183 | | sage: K.find(1) |
184 | | 1 |
185 | | |
186 | | Now we put them in the same cell:: |
187 | | |
188 | | sage: K.union_find(0,1) |
189 | | sage: K.find(0) |
190 | | 0 |
191 | | sage: K.find(1) |
192 | | 0 |
193 | | """ |
194 | | self._union_find(a, b) |
195 | | |
196 | | cdef int _union_find(self, int a, int b): |
197 | | cdef int aRoot, bRoot |
198 | | aRoot = self._find(a) |
199 | | bRoot = self._find(b) |
200 | | self._union_roots(aRoot, bRoot) |
201 | | return aRoot != bRoot |
202 | | |
203 | | cdef void _union_roots(self, int a, int b): |
204 | | if a < b: |
205 | | self.elements[b] = a |
206 | | self.sizes[b] += self.sizes[a] |
207 | | elif a > b: |
208 | | self.elements[a] = b |
209 | | self.sizes[a] += self.sizes[b] |
210 | | |
211 | | def is_finer_than(self, other, n): |
212 | | """ |
213 | | Partition P is finer than partition Q if every cell of P is a |
214 | | subset of a cell of Q. |
215 | | |
216 | | EXAMPLE:: |
217 | | |
218 | | sage: from sage.graphs.graph_isom import OrbitPartition |
219 | | sage: K = OrbitPartition(20) |
220 | | sage: K.find(7) |
221 | | 7 |
222 | | sage: K.union_find(7, 12) |
223 | | sage: K.find(12) |
224 | | 7 |
225 | | sage: J = OrbitPartition(20) |
226 | | sage: J.is_finer_than(K, 20) |
227 | | True |
228 | | sage: K.is_finer_than(J, 20) |
229 | | False |
230 | | """ |
231 | | return self._is_finer_than(other, n) == 1 |
232 | | |
233 | | cdef int _is_finer_than(self, OrbitPartition other, int n): |
234 | | cdef int i |
235 | | for i from 0 <= i < n: |
236 | | if self.elements[i] != -1 and other.find(self.find(i)) != other.find(i): |
237 | | return 0 |
238 | | return 1 |
239 | | |
240 | | def vee_with(self, other, n): |
241 | | """ |
242 | | Merges the minimal number of cells such that other is finer than |
243 | | self. |
244 | | |
245 | | EXAMPLE:: |
246 | | |
247 | | sage: from sage.graphs.graph_isom import OrbitPartition |
248 | | sage: K = OrbitPartition(20) |
249 | | sage: K.union_find(7, 12) |
250 | | sage: J = OrbitPartition(20) |
251 | | sage: J.is_finer_than(K, 20) |
252 | | True |
253 | | sage: K.is_finer_than(J, 20) |
254 | | False |
255 | | sage: J.vee_with(K, 20) |
256 | | sage: K.is_finer_than(J, 20) |
257 | | True |
258 | | """ |
259 | | self._vee_with(other, n) |
260 | | |
261 | | cdef void _vee_with(self, OrbitPartition other, int n): |
262 | | cdef int i |
263 | | for i from 0 <= i < n: |
264 | | if self.elements[i] == -1: |
265 | | self._union_roots(i, self.find(other.find(i))) |
266 | | |
267 | | cdef int _is_min_cell_rep(self, int i): |
268 | | if self.elements[i] == -1: |
269 | | return 1 |
270 | | return 0 |
271 | | |
272 | | cdef int _is_fixed(self, int i): |
273 | | if self.elements[i] == -1 and self.sizes[i] == 1: |
274 | | return 1 |
275 | | return 0 |
276 | | |
277 | | def incorporate_permutation(self, gamma): |
278 | | """ |
279 | | Unions the cells of self which contain common elements of some |
280 | | orbit of gamma. |
281 | | |
282 | | INPUT: |
283 | | |
284 | | |
285 | | - ``gamma`` - a permutation, in list notation |
286 | | |
287 | | |
288 | | EXAMPLE:: |
289 | | |
290 | | sage: from sage.graphs.graph_isom import OrbitPartition |
291 | | sage: O = OrbitPartition(9) |
292 | | sage: O.incorporate_permutation([0,1,3,2,5,6,7,4,8]) |
293 | | sage: for i in range(9): |
294 | | ... print i, O.find(i) |
295 | | 0 0 |
296 | | 1 1 |
297 | | 2 2 |
298 | | 3 2 |
299 | | 4 4 |
300 | | 5 4 |
301 | | 6 4 |
302 | | 7 4 |
303 | | 8 8 |
304 | | """ |
305 | | cdef int k, n = len(gamma) |
306 | | cdef int *_gamma = <int *> sage_malloc( n * sizeof(int) ) |
307 | | if not _gamma: |
308 | | raise MemoryError("Error allocating memory.") |
309 | | for k from 0 <= k < n: |
310 | | _gamma[k] = gamma[k] |
311 | | self._incorporate_permutation(_gamma, n) |
312 | | sage_free(_gamma) |
313 | | |
314 | | cdef int _incorporate_permutation(self, int *gamma, int n): |
315 | | cdef int i, ch = 0, k |
316 | | for i from 0 <= i < n: |
317 | | k = self._union_find(i, gamma[i]) |
318 | | if (not ch) and k: |
319 | | ch = 1 |
320 | | return ch |
321 | | |
322 | | cdef OrbitPartition orbit_partition_from_list_perm(int *gamma, int n): |
323 | | cdef int i |
324 | | cdef OrbitPartition O |
325 | | O = OrbitPartition(n) |
326 | | for i from 0 <= i < n: |
327 | | if i != gamma[i]: |
328 | | O._union_find(i, gamma[i]) |
329 | | return O |
330 | | |
331 | | cdef class PartitionStack: |
332 | | """ |
333 | | TODO: documentation |
334 | | |
335 | | EXAMPLES:: |
336 | | |
337 | | sage: from sage.graphs.graph_isom import PartitionStack |
338 | | sage: from sage.graphs.base.sparse_graph import SparseGraph |
339 | | sage: P = PartitionStack([range(9, -1, -1)]) |
340 | | sage: P.set_k(1) |
341 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
342 | | 0 |
343 | | sage: P.set_k(2) |
344 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
345 | | 0 |
346 | | sage: P.set_k(3) |
347 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
348 | | 4 |
349 | | sage: P.set_k(4) |
350 | | sage: P.sort_by_function(0, [0,1], 10) |
351 | | 0 |
352 | | sage: P.set_k(5) |
353 | | sage: P.sort_by_function(2, [1,0], 10) |
354 | | 2 |
355 | | sage: P.set_k(6) |
356 | | sage: P.sort_by_function(4, [1,0], 10) |
357 | | 4 |
358 | | sage: P.set_k(7) |
359 | | sage: P.sort_by_function(6, [1,0], 10) |
360 | | 6 |
361 | | sage: P |
362 | | (5,9,7,1,6,2,8,0,4,3) |
363 | | (5,9,7,1|6,2,8,0|4|3) |
364 | | (5,9|7,1|6,2,8,0|4|3) |
365 | | (5,9|7,1|6,2|8,0|4|3) |
366 | | (5|9|7,1|6,2|8,0|4|3) |
367 | | (5|9|7|1|6,2|8,0|4|3) |
368 | | (5|9|7|1|6|2|8,0|4|3) |
369 | | (5|9|7|1|6|2|8|0|4|3) |
370 | | sage: P.is_discrete() |
371 | | 1 |
372 | | sage: P.set_k(6) |
373 | | sage: P.is_discrete() |
374 | | 0 |
375 | | |
376 | | :: |
377 | | |
378 | | sage: G = SparseGraph(10) |
379 | | sage: for i,j,_ in graphs.PetersenGraph().edge_iterator(): |
380 | | ... G.add_arc(i,j) |
381 | | ... G.add_arc(j,i) |
382 | | sage: P = PartitionStack(10) |
383 | | sage: P.set_k(1) |
384 | | sage: P.split_vertex(0) |
385 | | sage: P.refine(G, [0], 10, 0, 1) |
386 | | sage: P |
387 | | (0,2,3,6,7,8,9,1,4,5) |
388 | | (0|2,3,6,7,8,9|1,4,5) |
389 | | sage: P.set_k(2) |
390 | | sage: P.split_vertex(1) |
391 | | sage: P.refine(G, [7], 10, 0, 1) |
392 | | sage: P |
393 | | (0,3,7,8,9,2,6,1,4,5) |
394 | | (0|3,7,8,9,2,6|1,4,5) |
395 | | (0|3,7,8,9|2,6|1|4,5) |
396 | | """ |
397 | | def __new__(self, data): |
398 | | cdef int j, k, n |
399 | | cdef PartitionStack _data |
400 | | try: |
401 | | n = int(data) |
402 | | self.entries = <int *> sage_malloc( n * sizeof(int) ) |
403 | | if not self.entries: |
404 | | raise MemoryError("Error allocating memory.") |
405 | | self.levels = <int *> sage_malloc( n * sizeof(int) ) |
406 | | if not self.levels: |
407 | | sage_free(self.entries) |
408 | | raise MemoryError("Error allocating memory.") |
409 | | for k from 0 <= k < n-1: |
410 | | self.entries[k] = k |
411 | | self.levels[k] = n |
412 | | self.entries[n-1] = n-1 |
413 | | self.levels[n-1] = -1 |
414 | | self.k = 0 |
415 | | except: |
416 | | if isinstance(data, list): |
417 | | n = sum([len(datum) for datum in data]) |
418 | | self.entries = <int *> sage_malloc( n * sizeof(int) ) |
419 | | if not self.entries: |
420 | | raise MemoryError("Error allocating memory.") |
421 | | self.levels = <int *> sage_malloc( n * sizeof(int) ) |
422 | | if not self.levels: |
423 | | sage_free(self.entries) |
424 | | raise MemoryError("Error allocating memory.") |
425 | | j = 0 |
426 | | k = 0 |
427 | | for cell in data: |
428 | | for entry in cell: |
429 | | self.entries[j] = entry |
430 | | self.levels[j] = n |
431 | | j += 1 |
432 | | self.levels[j-1] = 0 |
433 | | self._percolate(k, j-1) |
434 | | k = j |
435 | | self.levels[j-1] = -1 |
436 | | self.k = 0 |
437 | | elif isinstance(data, PartitionStack): |
438 | | _data = data |
439 | | j = 0 |
440 | | while _data.levels[j] != -1: j += 1 |
441 | | n = j + 1 |
442 | | self.entries = <int *> sage_malloc( n * sizeof(int) ) |
443 | | if not self.entries: |
444 | | raise MemoryError("Error allocating memory.") |
445 | | self.levels = <int *> sage_malloc( n * sizeof(int) ) |
446 | | if not self.levels: |
447 | | sage_free(self.entries) |
448 | | raise MemoryError("Error allocating memory.") |
449 | | for k from 0 <= k < n: |
450 | | self.entries[k] = _data.entries[k] |
451 | | self.levels[k] = _data.levels[k] |
452 | | self.k = _data.k |
453 | | else: |
454 | | raise ValueError("Input must be an int, a list of lists, or a PartitionStack.") |
455 | | |
456 | | def __dealloc__(self): |
457 | | sage_free(self.entries) |
458 | | sage_free(self.levels) |
459 | | |
460 | | def __repr__(self): |
461 | | """ |
462 | | EXAMPLE:: |
463 | | |
464 | | sage: from sage.graphs.graph_isom import PartitionStack |
465 | | sage: P = PartitionStack([range(9, -1, -1)]) |
466 | | sage: P.set_k(1) |
467 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
468 | | 0 |
469 | | sage: P.set_k(2) |
470 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
471 | | 0 |
472 | | sage: P.set_k(3) |
473 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
474 | | 4 |
475 | | sage: P.set_k(4) |
476 | | sage: P.sort_by_function(0, [0,1], 10) |
477 | | 0 |
478 | | sage: P.set_k(5) |
479 | | sage: P.sort_by_function(2, [1,0], 10) |
480 | | 2 |
481 | | sage: P.set_k(6) |
482 | | sage: P.sort_by_function(4, [1,0], 10) |
483 | | 4 |
484 | | sage: P.set_k(7) |
485 | | sage: P.sort_by_function(6, [1,0], 10) |
486 | | 6 |
487 | | sage: P |
488 | | (5,9,7,1,6,2,8,0,4,3) |
489 | | (5,9,7,1|6,2,8,0|4|3) |
490 | | (5,9|7,1|6,2,8,0|4|3) |
491 | | (5,9|7,1|6,2|8,0|4|3) |
492 | | (5|9|7,1|6,2|8,0|4|3) |
493 | | (5|9|7|1|6,2|8,0|4|3) |
494 | | (5|9|7|1|6|2|8,0|4|3) |
495 | | (5|9|7|1|6|2|8|0|4|3) |
496 | | """ |
497 | | k = 0 |
498 | | s = '' |
499 | | while (k == 0 or self.levels[k-1] != -1) and k <= self.k: |
500 | | s += self.repr_at_k(k) + '\n' |
501 | | k += 1 |
502 | | return s |
503 | | |
504 | | def repr_at_k(self, k): |
505 | | """ |
506 | | Return the k-th line of the representation of self, i.e. the k-th |
507 | | partition in the stack. |
508 | | |
509 | | EXAMPLE:: |
510 | | |
511 | | sage: from sage.graphs.graph_isom import PartitionStack |
512 | | sage: P = PartitionStack([range(9, -1, -1)]) |
513 | | sage: P.set_k(1) |
514 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
515 | | 0 |
516 | | sage: P.set_k(2) |
517 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
518 | | 0 |
519 | | sage: P.set_k(3) |
520 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
521 | | 4 |
522 | | sage: P.set_k(4) |
523 | | sage: P.sort_by_function(0, [0,1], 10) |
524 | | 0 |
525 | | sage: P.set_k(5) |
526 | | sage: P.sort_by_function(2, [1,0], 10) |
527 | | 2 |
528 | | sage: P.set_k(6) |
529 | | sage: P.sort_by_function(4, [1,0], 10) |
530 | | 4 |
531 | | sage: P.set_k(7) |
532 | | sage: P.sort_by_function(6, [1,0], 10) |
533 | | 6 |
534 | | sage: P |
535 | | (5,9,7,1,6,2,8,0,4,3) |
536 | | (5,9,7,1|6,2,8,0|4|3) |
537 | | (5,9|7,1|6,2,8,0|4|3) |
538 | | (5,9|7,1|6,2|8,0|4|3) |
539 | | (5|9|7,1|6,2|8,0|4|3) |
540 | | (5|9|7|1|6,2|8,0|4|3) |
541 | | (5|9|7|1|6|2|8,0|4|3) |
542 | | (5|9|7|1|6|2|8|0|4|3) |
543 | | |
544 | | :: |
545 | | |
546 | | sage: P.repr_at_k(0) |
547 | | '(5,9,7,1,6,2,8,0,4,3)' |
548 | | sage: P.repr_at_k(1) |
549 | | '(5,9,7,1|6,2,8,0|4|3)' |
550 | | """ |
551 | | s = '(' |
552 | | i = 0 |
553 | | while i == 0 or self.levels[i-1] != -1: |
554 | | s += str(self.entries[i]) |
555 | | if self.levels[i] <= k: |
556 | | s += '|' |
557 | | else: |
558 | | s += ',' |
559 | | i += 1 |
560 | | s = s[:-1] + ')' |
561 | | return s |
562 | | |
563 | | def set_k(self, k): |
564 | | """ |
565 | | Sets self.k, the index of the finest partition. |
566 | | |
567 | | :: |
568 | | |
569 | | sage: from sage.graphs.graph_isom import PartitionStack |
570 | | sage: P = PartitionStack([range(9, -1, -1)]) |
571 | | sage: P.set_k(1) |
572 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
573 | | 0 |
574 | | sage: P.set_k(2) |
575 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
576 | | 0 |
577 | | sage: P.set_k(3) |
578 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
579 | | 4 |
580 | | sage: P.set_k(4) |
581 | | sage: P.sort_by_function(0, [0,1], 10) |
582 | | 0 |
583 | | sage: P.set_k(5) |
584 | | sage: P.sort_by_function(2, [1,0], 10) |
585 | | 2 |
586 | | sage: P.set_k(6) |
587 | | sage: P.sort_by_function(4, [1,0], 10) |
588 | | 4 |
589 | | sage: P.set_k(7) |
590 | | sage: P.sort_by_function(6, [1,0], 10) |
591 | | 6 |
592 | | sage: P |
593 | | (5,9,7,1,6,2,8,0,4,3) |
594 | | (5,9,7,1|6,2,8,0|4|3) |
595 | | (5,9|7,1|6,2,8,0|4|3) |
596 | | (5,9|7,1|6,2|8,0|4|3) |
597 | | (5|9|7,1|6,2|8,0|4|3) |
598 | | (5|9|7|1|6,2|8,0|4|3) |
599 | | (5|9|7|1|6|2|8,0|4|3) |
600 | | (5|9|7|1|6|2|8|0|4|3) |
601 | | |
602 | | :: |
603 | | |
604 | | sage: P.set_k(2) |
605 | | sage: P |
606 | | (5,9,7,1,6,2,8,0,4,3) |
607 | | (5,9,7,1|6,2,8,0|4|3) |
608 | | (5,9|7,1|6,2,8,0|4|3) |
609 | | """ |
610 | | self.k = k |
611 | | |
612 | | def is_discrete(self): |
613 | | """ |
614 | | Returns whether the partition consists of only singletons. |
615 | | |
616 | | EXAMPLE:: |
617 | | |
618 | | sage: from sage.graphs.graph_isom import PartitionStack |
619 | | sage: P = PartitionStack([range(9, -1, -1)]) |
620 | | sage: P.set_k(1) |
621 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
622 | | 0 |
623 | | sage: P.set_k(2) |
624 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
625 | | 0 |
626 | | sage: P.set_k(3) |
627 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
628 | | 4 |
629 | | sage: P.set_k(4) |
630 | | sage: P.sort_by_function(0, [0,1], 10) |
631 | | 0 |
632 | | sage: P.set_k(5) |
633 | | sage: P.sort_by_function(2, [1,0], 10) |
634 | | 2 |
635 | | sage: P.set_k(6) |
636 | | sage: P.sort_by_function(4, [1,0], 10) |
637 | | 4 |
638 | | sage: P.set_k(7) |
639 | | sage: P.sort_by_function(6, [1,0], 10) |
640 | | 6 |
641 | | sage: P |
642 | | (5,9,7,1,6,2,8,0,4,3) |
643 | | (5,9,7,1|6,2,8,0|4|3) |
644 | | (5,9|7,1|6,2,8,0|4|3) |
645 | | (5,9|7,1|6,2|8,0|4|3) |
646 | | (5|9|7,1|6,2|8,0|4|3) |
647 | | (5|9|7|1|6,2|8,0|4|3) |
648 | | (5|9|7|1|6|2|8,0|4|3) |
649 | | (5|9|7|1|6|2|8|0|4|3) |
650 | | sage: P.is_discrete() |
651 | | True |
652 | | sage: P.set_k(2) |
653 | | sage: P |
654 | | (5,9,7,1,6,2,8,0,4,3) |
655 | | (5,9,7,1|6,2,8,0|4|3) |
656 | | (5,9|7,1|6,2,8,0|4|3) |
657 | | sage: P.is_discrete() |
658 | | False |
659 | | """ |
660 | | return (self._is_discrete() == 1) |
661 | | |
662 | | cdef int _is_discrete(self): |
663 | | cdef int i = 0 |
664 | | while True: |
665 | | if self.levels[i] > self.k: |
666 | | return 0 |
667 | | if self.levels[i] == -1: break |
668 | | i += 1 |
669 | | return 1 |
670 | | |
671 | | def num_cells(self): |
672 | | """ |
673 | | Return the number of cells in the finest partition. |
674 | | |
675 | | EXAMPLE:: |
676 | | |
677 | | sage: from sage.graphs.graph_isom import PartitionStack |
678 | | sage: P = PartitionStack([range(9, -1, -1)]) |
679 | | sage: P.set_k(1) |
680 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
681 | | 0 |
682 | | sage: P |
683 | | (1,9,7,5,0,2,8,6,4,3) |
684 | | (1,9,7,5|0,2,8,6|4|3) |
685 | | sage: P.num_cells() |
686 | | 4 |
687 | | sage: P.set_k(2) |
688 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
689 | | 0 |
690 | | sage: P |
691 | | (5,9,1,7,0,2,8,6,4,3) |
692 | | (5,9,1,7|0,2,8,6|4|3) |
693 | | (5,9|1,7|0,2,8,6|4|3) |
694 | | sage: P.num_cells() |
695 | | 5 |
696 | | """ |
697 | | return self._num_cells() |
698 | | |
699 | | cdef int _num_cells(self): |
700 | | cdef int i = 0, j = 1 |
701 | | while self.levels[i] != -1: |
702 | | #for i from 0 <= i < n-1: |
703 | | if self.levels[i] <= self.k: |
704 | | j += 1 |
705 | | i += 1 |
706 | | return j |
707 | | |
708 | | def sat_225(self, n): |
709 | | """ |
710 | | Whether the finest partition satisfies the hypotheses of Lemma 2.25 |
711 | | in [1]. |
712 | | |
713 | | EXAMPLE:: |
714 | | |
715 | | sage: from sage.graphs.graph_isom import PartitionStack |
716 | | sage: P = PartitionStack([range(9, -1, -1)]) |
717 | | sage: P |
718 | | (0,9,8,7,6,5,4,3,2,1) |
719 | | sage: P.sat_225(10) |
720 | | False |
721 | | sage: P.set_k(1) |
722 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
723 | | 0 |
724 | | sage: P |
725 | | (1,9,7,5,0,2,8,6,4,3) |
726 | | (1,9,7,5|0,2,8,6|4|3) |
727 | | sage: P.sat_225(10) |
728 | | False |
729 | | sage: P.set_k(2) |
730 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
731 | | 0 |
732 | | sage: P |
733 | | (5,9,1,7,0,2,8,6,4,3) |
734 | | (5,9,1,7|0,2,8,6|4|3) |
735 | | (5,9|1,7|0,2,8,6|4|3) |
736 | | sage: P.sat_225(10) |
737 | | False |
738 | | sage: P.set_k(3) |
739 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
740 | | 4 |
741 | | sage: P |
742 | | (5,9,1,7,2,6,0,8,4,3) |
743 | | (5,9,1,7|2,6,0,8|4|3) |
744 | | (5,9|1,7|2,6,0,8|4|3) |
745 | | (5,9|1,7|2,6|0,8|4|3) |
746 | | sage: P.sat_225(10) |
747 | | True |
748 | | """ |
749 | | return self._sat_225(n) == 1 |
750 | | |
751 | | cdef int _sat_225(self, int n): |
752 | | cdef int i, in_cell = 0 |
753 | | cdef int nontrivial_cells = 0 |
754 | | cdef int total_cells = self._num_cells() |
755 | | if n <= total_cells + 4: |
756 | | return 1 |
757 | | for i from 0 <= i < n-1: |
758 | | if self.levels[i] <= self.k: |
759 | | if in_cell: |
760 | | nontrivial_cells += 1 |
761 | | in_cell = 0 |
762 | | else: |
763 | | in_cell = 1 |
764 | | if in_cell: |
765 | | nontrivial_cells += 1 |
766 | | if n == total_cells + nontrivial_cells: |
767 | | return 1 |
768 | | if n == total_cells + nontrivial_cells + 1: |
769 | | return 1 |
770 | | return 0 |
771 | | |
772 | | cdef int _is_min_cell_rep(self, int i, int k): |
773 | | return i == 0 or self.levels[i-1] <= k |
774 | | |
775 | | cdef int _is_fixed(self, int i, int k): |
776 | | """ |
777 | | Assuming you already know it is a minimum cell representative. |
778 | | """ |
779 | | return self.levels[i] <= k |
780 | | |
781 | | def split_vertex(self, v): |
782 | | """ |
783 | | Splits the cell in self(k) containing v, putting new cells in place |
784 | | in self(k). |
785 | | |
786 | | EXAMPLE:: |
787 | | |
788 | | sage: from sage.graphs.graph_isom import PartitionStack |
789 | | sage: P = PartitionStack([range(9, -1, -1)]) |
790 | | sage: P |
791 | | (0,9,8,7,6,5,4,3,2,1) |
792 | | sage: P.split_vertex(2) |
793 | | sage: P |
794 | | (2|0,9,8,7,6,5,4,3,1) |
795 | | """ |
796 | | self._split_vertex(v) |
797 | | |
798 | | cdef int _split_vertex(self, int v): |
799 | | cdef int i = 0, j |
800 | | while self.entries[i] != v: |
801 | | i += 1 |
802 | | j = i |
803 | | while self.levels[i] > self.k: |
804 | | i += 1 |
805 | | if j == 0 or self.levels[j-1] <= self.k: |
806 | | self._percolate(j+1, i) |
807 | | else: |
808 | | while j != 0 and self.levels[j-1] > self.k: |
809 | | self.entries[j] = self.entries[j-1] |
810 | | j -= 1 |
811 | | self.entries[j] = v |
812 | | self.levels[j] = self.k |
813 | | return j |
814 | | |
815 | | def percolate(self, start, end): |
816 | | """ |
817 | | Perform one round of bubble sort, moving the smallest element to |
818 | | the front. |
819 | | |
820 | | EXAMPLE:: |
821 | | |
822 | | sage: from sage.graphs.graph_isom import PartitionStack |
823 | | sage: P = PartitionStack([range(9, -1, -1)]) |
824 | | sage: P |
825 | | (0,9,8,7,6,5,4,3,2,1) |
826 | | sage: P.percolate(2,7) |
827 | | sage: P |
828 | | (0,9,3,8,7,6,5,4,2,1) |
829 | | """ |
830 | | self._percolate(start, end) |
831 | | |
832 | | cdef void _percolate(self, int start, int end): |
833 | | cdef int i, temp |
834 | | for i from end >= i > start: |
835 | | if self.entries[i] < self.entries[i-1]: |
836 | | temp = self.entries[i] |
837 | | self.entries[i] = self.entries[i-1] |
838 | | self.entries[i-1] = temp |
839 | | |
840 | | def sort_by_function(self, start, degrees, n): |
841 | | """ |
842 | | Sort the cell starting at start using a counting sort, where |
843 | | degrees is the function giving the sort. Result is the cell is |
844 | | subdivided into cells which have elements all of the same 'degree,' |
845 | | in order. |
846 | | |
847 | | EXAMPLE:: |
848 | | |
849 | | sage: from sage.graphs.graph_isom import PartitionStack |
850 | | sage: P = PartitionStack([range(9, -1, -1)]) |
851 | | sage: P.set_k(1) |
852 | | sage: P |
853 | | (0,9,8,7,6,5,4,3,2,1) |
854 | | (0,9,8,7,6,5,4,3,2,1) |
855 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
856 | | 0 |
857 | | sage: P |
858 | | (1,9,7,5,0,2,8,6,4,3) |
859 | | (1,9,7,5|0,2,8,6|4|3) |
860 | | """ |
861 | | cdef int i |
862 | | cdef int *degs = <int *> sage_malloc( ( 3 * n + 1 ) * sizeof(int) ) |
863 | | if not degs: |
864 | | raise MemoryError("Couldn't allocate...") |
865 | | for i from 0 <= i < len(degrees): |
866 | | degs[i] = degrees[i] |
867 | | X = self._sort_by_function(start, degs, n) |
868 | | sage_free(degs) |
869 | | return X |
870 | | |
871 | | cdef int _sort_by_function(self, int start, int *degrees, int n): |
872 | | cdef int i, j, m = 2*n, max, max_location |
873 | | cdef int *counts = degrees + n, *output = degrees + 2*n + 1 |
874 | | # print '|'.join(['%02d'%self.entries[iii] for iii in range(n)]) |
875 | | # print '|'.join(['%02d'%self.levels[iii] for iii in range(n)]) |
876 | | # print '|'.join(['%02d'%degrees[iii] for iii in range(n)]) |
877 | | # print '|'.join(['%02d'%counts[iii] for iii in range(n)]) |
878 | | # print '|'.join(['%02d'%output[iii] for iii in range(n)]) |
879 | | |
880 | | for i from 0 <= i <= n: |
881 | | counts[i] = 0 |
882 | | i = 0 |
883 | | while self.levels[i+start] > self.k: |
884 | | counts[degrees[i]] += 1 |
885 | | i += 1 |
886 | | counts[degrees[i]] += 1 |
887 | | |
888 | | # i+start is the right endpoint of the cell now |
889 | | max = counts[0] |
890 | | max_location = 0 |
891 | | for j from 0 < j <= n: |
892 | | if counts[j] > max: |
893 | | max = counts[j] |
894 | | max_location = j |
895 | | counts[j] += counts[j - 1] |
896 | | |
897 | | for j from i >= j >= 0: |
898 | | counts[degrees[j]] -= 1 |
899 | | output[counts[degrees[j]]] = self.entries[start+j] |
900 | | |
901 | | max_location = counts[max_location]+start |
902 | | |
903 | | for j from 0 <= j <= i: |
904 | | self.entries[start+j] = output[j] |
905 | | |
906 | | j = 1 |
907 | | while j <= n and counts[j] <= i: |
908 | | if counts[j] > 0: |
909 | | self.levels[start + counts[j] - 1] = self.k |
910 | | self._percolate(start + counts[j-1], start + counts[j] - 1) |
911 | | j += 1 |
912 | | |
913 | | return max_location |
914 | | |
915 | | def clear(self): |
916 | | """ |
917 | | Merges all cells in the partition stack. |
918 | | |
919 | | EXAMPLE:: |
920 | | |
921 | | sage: from sage.graphs.graph_isom import PartitionStack |
922 | | sage: P = PartitionStack([range(9, -1, -1)]) |
923 | | sage: P.set_k(1) |
924 | | sage: P |
925 | | (0,9,8,7,6,5,4,3,2,1) |
926 | | (0,9,8,7,6,5,4,3,2,1) |
927 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
928 | | 0 |
929 | | sage: P |
930 | | (1,9,7,5,0,2,8,6,4,3) |
931 | | (1,9,7,5|0,2,8,6|4|3) |
932 | | sage: P |
933 | | (1,9,7,5,0,2,8,6,4,3) |
934 | | (1,9,7,5|0,2,8,6|4|3) |
935 | | sage: P.clear() |
936 | | sage: P |
937 | | (1,9,7,5,0,2,8,6,4,3) |
938 | | (1,9,7,5,0,2,8,6,4,3) |
939 | | """ |
940 | | self._clear() |
941 | | |
942 | | cdef void _clear(self): |
943 | | cdef int i = 0, j = 0 |
944 | | while self.levels[i] != -1: |
945 | | if self.levels[i] >= self.k: |
946 | | self.levels[i] += 1 |
947 | | if self.levels[i] < self.k: |
948 | | self._percolate(j, i) |
949 | | j = i + 1 |
950 | | i+=1 |
951 | | |
952 | | def refine(self, CGraph G, alpha, n, dig, uif, test=False): |
953 | | """ |
954 | | Implementation of Algorithm 2.5 in [1]. |
955 | | |
956 | | EXAMPLE:: |
957 | | |
958 | | sage: from sage.graphs.graph_isom import PartitionStack |
959 | | sage: from sage.graphs.base.sparse_graph import SparseGraph |
960 | | sage: P = PartitionStack([range(9, -1, -1)]) |
961 | | sage: P.set_k(1) |
962 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
963 | | 0 |
964 | | sage: P.set_k(2) |
965 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
966 | | 0 |
967 | | sage: P.set_k(3) |
968 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
969 | | 4 |
970 | | sage: P.set_k(4) |
971 | | sage: P.sort_by_function(0, [0,1], 10) |
972 | | 0 |
973 | | sage: P.set_k(5) |
974 | | sage: P.sort_by_function(2, [1,0], 10) |
975 | | 2 |
976 | | sage: P.set_k(6) |
977 | | sage: P.sort_by_function(4, [1,0], 10) |
978 | | 4 |
979 | | sage: P.set_k(7) |
980 | | sage: P.sort_by_function(6, [1,0], 10) |
981 | | 6 |
982 | | sage: P |
983 | | (5,9,7,1,6,2,8,0,4,3) |
984 | | (5,9,7,1|6,2,8,0|4|3) |
985 | | (5,9|7,1|6,2,8,0|4|3) |
986 | | (5,9|7,1|6,2|8,0|4|3) |
987 | | (5|9|7,1|6,2|8,0|4|3) |
988 | | (5|9|7|1|6,2|8,0|4|3) |
989 | | (5|9|7|1|6|2|8,0|4|3) |
990 | | (5|9|7|1|6|2|8|0|4|3) |
991 | | sage: P.is_discrete() |
992 | | 1 |
993 | | sage: P.set_k(6) |
994 | | sage: P.is_discrete() |
995 | | 0 |
996 | | |
997 | | :: |
998 | | |
999 | | sage: G = SparseGraph(10) |
1000 | | sage: for i,j,_ in graphs.PetersenGraph().edge_iterator(): |
1001 | | ... G.add_arc(i,j) |
1002 | | ... G.add_arc(j,i) |
1003 | | sage: P = PartitionStack(10) |
1004 | | sage: P.set_k(1) |
1005 | | sage: P.split_vertex(0) |
1006 | | sage: P.refine(G, [0], 10, 0, 1) |
1007 | | sage: P |
1008 | | (0,2,3,6,7,8,9,1,4,5) |
1009 | | (0|2,3,6,7,8,9|1,4,5) |
1010 | | sage: P.set_k(2) |
1011 | | sage: P.split_vertex(1) |
1012 | | sage: P.refine(G, [7], 10, 0, 1) |
1013 | | sage: P |
1014 | | (0,3,7,8,9,2,6,1,4,5) |
1015 | | (0|3,7,8,9,2,6|1,4,5) |
1016 | | (0|3,7,8,9|2,6|1|4,5) |
1017 | | """ |
1018 | | cdef int *_alpha, i, j |
1019 | | _alpha = <int *> sage_malloc( ( 4 * n + 1 )* sizeof(int) ) |
1020 | | if not _alpha: |
1021 | | raise MemoryError("Memory!") |
1022 | | for i from 0 <= i < len(alpha): |
1023 | | _alpha[i] = alpha[i] |
1024 | | _alpha[len(alpha)] = -1 |
1025 | | if test: |
1026 | | self.test_refine(_alpha, n, G, dig, uif) |
1027 | | else: |
1028 | | self._refine(_alpha, n, G, dig, uif) |
1029 | | sage_free(_alpha) |
1030 | | |
1031 | | cdef int test_refine(self, int *alpha, int n, CGraph g, int dig, int uif) except? -1: |
1032 | | cdef int i, j, result |
1033 | | initial_partition = [] # this includes the vertex just split out... |
1034 | | i = 0 |
1035 | | cell = [] |
1036 | | while i < n: |
1037 | | cell.append(self.entries[i]) |
1038 | | while self.levels[i] > self.k: |
1039 | | i += 1 |
1040 | | cell.append(self.entries[i]) |
1041 | | i += 1 |
1042 | | initial_partition.append(cell) |
1043 | | cell = [] |
1044 | | # |
1045 | | result = self._refine(alpha, n, g, dig, uif) |
1046 | | # |
1047 | | terminal_partition = [] |
1048 | | i = 0 |
1049 | | cell = [] |
1050 | | while i < n: |
1051 | | cell.append(self.entries[i]) |
1052 | | while self.levels[i] > self.k: |
1053 | | i += 1 |
1054 | | cell.append(self.entries[i]) |
1055 | | i += 1 |
1056 | | terminal_partition.append(cell) |
1057 | | cell = [] |
1058 | | # |
1059 | | if dig: |
1060 | | G = DiGraph(n, loops=True) |
1061 | | else: |
1062 | | G = Graph(n) |
1063 | | for i from 0 <= i < n: |
1064 | | for j from 0 <= j < n: |
1065 | | if g.has_arc_unsafe(i, j): |
1066 | | G.add_edge(i,j) |
1067 | | verify_partition_refinement(G, initial_partition, terminal_partition) |
1068 | | return result |
1069 | | |
1070 | | cdef int _refine(self, int *alpha, int n, CGraph G, int dig, int uif): |
1071 | | cdef int m = 0, j # - m iterates through alpha, the indicator cells |
1072 | | # - j iterates through the cells of the partition |
1073 | | cdef int i, t, s, r # local variables: |
1074 | | # - s plays a double role: outer role indicates whether |
1075 | | # splitting the cell is necessary, inner role is as an index |
1076 | | # for augmenting _alpha |
1077 | | # - i, r iterators |
1078 | | # - t: holds the first largest subcell from sort function |
1079 | | cdef int invariant = 1 |
1080 | | # as described in [1], an indicator function Lambda(G, Pi, nu) is |
1081 | | # needed to differentiate nonisomorphic branches on the search |
1082 | | # tree. The condition is simply that this invariant not depend |
1083 | | # on a simultaneous relabeling of the graph G, the root partition |
1084 | | # Pi, and the partition nest nu. Since the function will execute |
1085 | | # exactly the same way regardless of the labelling, anything that |
1086 | | # does not depend on self.entries goes... at least, anything cheap |
1087 | | cdef int *degrees = alpha + n # alpha assumed to be length 4*n + 1 for |
1088 | | # extra scratch space |
1089 | | while not self._is_discrete() and alpha[m] != -1: |
1090 | | invariant += 1 |
1091 | | j = 0 |
1092 | | while j < n: # j still points at a valid cell |
1093 | | invariant += 50 |
1094 | | # print ' ' |
1095 | | # print '|'.join(['%02d'%self.entries[iii] for iii in range(n)]) |
1096 | | # print '|'.join(['%02d'%self.levels[iii] for iii in range(n)]) |
1097 | | # print '|'.join(['%02d'%alpha[iii] for iii in range(n)]) |
1098 | | # print '|'.join(['%02d'%degrees[iii] for iii in range(n)]) |
1099 | | # print 'j =', j |
1100 | | # print 'm =', m |
1101 | | i = j; s = 0 |
1102 | | while True: |
1103 | | degrees[i-j] = self._degree(G, i, alpha[m]) |
1104 | | if degrees[i-j] != degrees[0]: s = 1 |
1105 | | i += 1 |
1106 | | if self.levels[i-1] <= self.k: break |
1107 | | # print '|'.join(['%02d'%degrees[iii] for iii in range(n)]) |
1108 | | # now: j points to this cell, |
1109 | | # i points to the next cell (before refinement) |
1110 | | if s: |
1111 | | invariant += 10 |
1112 | | t = self._sort_by_function(j, degrees, n) |
1113 | | # t now points to the first largest subcell |
1114 | | invariant += t |
1115 | | s = m |
1116 | | while alpha[s] != -1: |
1117 | | if alpha[s] == j: |
1118 | | alpha[s] = t |
1119 | | break |
1120 | | s += 1 |
1121 | | while alpha[s] != -1: s += 1 |
1122 | | r = j |
1123 | | while True: |
1124 | | if r == j or self.levels[r-1] == self.k: |
1125 | | if r != t: |
1126 | | alpha[s] = r |
1127 | | s += 1 |
1128 | | r += 1 |
1129 | | if r >= i: break |
1130 | | alpha[s] = -1 |
1131 | | invariant += self._degree(G, i-1, alpha[m]) |
1132 | | invariant += (i - j) |
1133 | | j = i |
1134 | | else: j = i |
1135 | | if not dig: m += 1; continue |
1136 | | # if we are looking at a digraph, also compute |
1137 | | # the reverse degrees and sort by them |
1138 | | j = 0 |
1139 | | while j < n: # j still points at a valid cell |
1140 | | invariant += 20 |
1141 | | # print ' ' |
1142 | | # print '|'.join(['%02d'%self.entries[iii] for iii in range(n)]) |
1143 | | # print '|'.join(['%02d'%self.levels[iii] for iii in range(n)]) |
1144 | | # print '|'.join(['%02d'%alpha[iii] for iii in range(n)]) |
1145 | | # print '|'.join(['%02d'%degrees[iii] for iii in range(n)]) |
1146 | | # print 'j =', j |
1147 | | # print 'm =', m |
1148 | | i = j; s = 0 |
1149 | | while True: |
1150 | | degrees[i-j] = self._degree_inv(G, i, alpha[m]) |
1151 | | if degrees[i-j] != degrees[0]: s = 1 |
1152 | | i += 1 |
1153 | | if self.levels[i-1] <= self.k: break |
1154 | | # now: j points to this cell, |
1155 | | # i points to the next cell (before refinement) |
1156 | | if s: |
1157 | | invariant += 7 |
1158 | | t = self._sort_by_function(j, degrees, n) |
1159 | | # t now points to the first largest subcell |
1160 | | invariant += t |
1161 | | s = m |
1162 | | while alpha[s] != -1: |
1163 | | if alpha[s] == j: |
1164 | | alpha[s] = t |
1165 | | break |
1166 | | s += 1 |
1167 | | while alpha[s] != -1: s += 1 |
1168 | | r = j |
1169 | | while True: |
1170 | | if r == j or self.levels[r-1] == self.k: |
1171 | | if r != t: |
1172 | | alpha[s] = r |
1173 | | s += 1 |
1174 | | r += 1 |
1175 | | if r >= i: break |
1176 | | alpha[s] = -1 |
1177 | | invariant += self._degree(G, i-1, alpha[m]) |
1178 | | invariant += (i - j) |
1179 | | j = i |
1180 | | else: j = i |
1181 | | m += 1 |
1182 | | if uif: |
1183 | | return invariant |
1184 | | else: |
1185 | | return 0 |
1186 | | |
1187 | | def degree(self, CGraph G, v, W): |
1188 | | """ |
1189 | | Returns the number of edges in G from self.entries[v] to a vertex |
1190 | | in W. |
1191 | | |
1192 | | EXAMPLE:: |
1193 | | |
1194 | | sage: from sage.graphs.graph_isom import PartitionStack |
1195 | | sage: from sage.graphs.base.sparse_graph import SparseGraph |
1196 | | sage: P = PartitionStack([range(9, -1, -1)]) |
1197 | | sage: P |
1198 | | (0,9,8,7,6,5,4,3,2,1) |
1199 | | sage: G = SparseGraph(10) |
1200 | | sage: G.add_arc(2,9) |
1201 | | sage: G.add_arc(3,9) |
1202 | | sage: G.add_arc(4,9) |
1203 | | sage: P.degree(G, 1, 0) |
1204 | | 3 |
1205 | | """ |
1206 | | cdef int j |
1207 | | j = self._degree(G, v, W) |
1208 | | return j |
1209 | | |
1210 | | cdef int _degree(self, CGraph G, int v, int W): |
1211 | | """ |
1212 | | G is a CGraph, and W points to the beginning of a cell in the k-th |
1213 | | part of the stack. |
1214 | | """ |
1215 | | cdef int i = 0 |
1216 | | v = self.entries[v] |
1217 | | while True: |
1218 | | if G.has_arc_unsafe(self.entries[W], v): |
1219 | | i += 1 |
1220 | | if self.levels[W] > self.k: W += 1 |
1221 | | else: break |
1222 | | return i |
1223 | | |
1224 | | cdef int _degree_inv(self, CGraph G, int v, int W): |
1225 | | """ |
1226 | | G is a CGraph, and W points to the beginning of a cell in the k-th |
1227 | | part of the stack. |
1228 | | """ |
1229 | | cdef int i = 0 |
1230 | | v = self.entries[v] |
1231 | | while True: |
1232 | | if G.has_arc_unsafe(v, self.entries[W]): |
1233 | | i += 1 |
1234 | | if self.levels[W] > self.k: W += 1 |
1235 | | else: break |
1236 | | return i |
1237 | | |
1238 | | cdef int _first_smallest_nontrivial(self, int *W, int n): |
1239 | | cdef int i = 0, j = 0, location = 0, min = n |
1240 | | while True: |
1241 | | W[i] = 0 |
1242 | | if self.levels[i] <= self.k: |
1243 | | if i != j and n > i - j + 1: |
1244 | | n = i - j + 1 |
1245 | | location = j |
1246 | | j = i + 1 |
1247 | | if self.levels[i] == -1: break |
1248 | | i += 1 |
1249 | | # location now points to the beginning of the first, smallest, |
1250 | | # nontrivial cell |
1251 | | while True: |
1252 | | if min > self.entries[location]: |
1253 | | min = self.entries[location] |
1254 | | W[self.entries[location]] = 1 |
1255 | | if self.levels[location] <= self.k: break |
1256 | | location += 1 |
1257 | | return min |
1258 | | |
1259 | | cdef void _get_permutation_from(self, PartitionStack zeta, int *gamma): |
1260 | | cdef int i = 0 |
1261 | | |
1262 | | while True: |
1263 | | gamma[zeta.entries[i]] = self.entries[i] |
1264 | | i += 1 |
1265 | | if self.levels[i-1] == -1: break |
1266 | | |
1267 | | cdef int _compare_with(self, CGraph G, int n, PartitionStack other): |
1268 | | cdef int i, j |
1269 | | for i from 0 <= i < n: |
1270 | | for j from 0 <= j < n: |
1271 | | if G.has_arc_unsafe(self.entries[i], self.entries[j]): |
1272 | | if not G.has_arc_unsafe(other.entries[i], other.entries[j]): |
1273 | | return 1 |
1274 | | elif G.has_arc_unsafe(other.entries[i], other.entries[j]): |
1275 | | return -1 |
1276 | | return 0 |
1277 | | |
1278 | | cdef int _is_automorphism(CGraph G, int n, int *gamma): |
1279 | | cdef int i, j |
1280 | | for i from 0 <= i < n: |
1281 | | for j from 0 <= j < n: |
1282 | | if G.has_arc_unsafe(i, j): |
1283 | | if not G.has_arc_unsafe(gamma[i], gamma[j]): |
1284 | | return 0 |
1285 | | return 1 |
1286 | | |
1287 | | def _term_pnest_graph(G, PartitionStack nu): |
1288 | | """ |
1289 | | BDM's G(nu): returns the graph G, relabeled in the order found in |
1290 | | nu[m], where m is the first index corresponding to a discrete |
1291 | | partition. Assumes nu is a terminal partition nest in T(G, Pi). |
1292 | | |
1293 | | EXAMPLE:: |
1294 | | |
1295 | | sage: from sage.graphs.graph_isom import PartitionStack |
1296 | | sage: from sage.graphs.base.sparse_graph import SparseGraph |
1297 | | sage: P = PartitionStack([range(9, -1, -1)]) |
1298 | | sage: P.set_k(1) |
1299 | | sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 10) |
1300 | | 0 |
1301 | | sage: P.set_k(2) |
1302 | | sage: P.sort_by_function(0, [2,1,2,1], 10) |
1303 | | 0 |
1304 | | sage: P.set_k(3) |
1305 | | sage: P.sort_by_function(4, [2,1,2,1], 10) |
1306 | | 4 |
1307 | | sage: P.set_k(4) |
1308 | | sage: P.sort_by_function(0, [0,1], 10) |
1309 | | 0 |
1310 | | sage: P.set_k(5) |
1311 | | sage: P.sort_by_function(2, [1,0], 10) |
1312 | | 2 |
1313 | | sage: P.set_k(6) |
1314 | | sage: P.sort_by_function(4, [1,0], 10) |
1315 | | 4 |
1316 | | sage: P.set_k(7) |
1317 | | sage: P.sort_by_function(6, [1,0], 10) |
1318 | | 6 |
1319 | | sage: P |
1320 | | (5,9,7,1,6,2,8,0,4,3) |
1321 | | (5,9,7,1|6,2,8,0|4|3) |
1322 | | (5,9|7,1|6,2,8,0|4|3) |
1323 | | (5,9|7,1|6,2|8,0|4|3) |
1324 | | (5|9|7,1|6,2|8,0|4|3) |
1325 | | (5|9|7|1|6,2|8,0|4|3) |
1326 | | (5|9|7|1|6|2|8,0|4|3) |
1327 | | (5|9|7|1|6|2|8|0|4|3) |
1328 | | sage: from sage.graphs.graph_isom import _term_pnest_graph |
1329 | | sage: _term_pnest_graph(graphs.PetersenGraph(), P).edges(labels=False) |
1330 | | [(0, 2), (0, 6), (0, 7), (1, 2), (1, 4), (1, 8), (2, 5), (3, 4), (3, 5), (3, 7), (4, 6), (5, 9), (6, 9), (7, 8), (8, 9)] |
1331 | | """ |
1332 | | cdef int i, j, n |
1333 | | cdef CGraph M |
1334 | | if isinstance(G, GenericGraph): |
1335 | | n = G.order() |
1336 | | H = G.copy() |
1337 | | else: # G is a CGraph |
1338 | | M = G |
1339 | | n = M.num_verts |
1340 | | if isinstance(G, SparseGraph): |
1341 | | H = SparseGraph(n) |
1342 | | else: |
1343 | | H = DenseGraph(n) |
1344 | | d = {} |
1345 | | for i from 0 <= i < n: |
1346 | | d[nu.entries[i]] = i |
1347 | | if isinstance(G, GenericGraph): |
1348 | | H.relabel(d) |
1349 | | else: |
1350 | | for i from 0 <= i < n: |
1351 | | for j in G.out_neighbors(i): |
1352 | | H.add_arc(d[i],d[j]) |
1353 | | return H |
1354 | | |
1355 | | def search_tree(G, Pi, lab=True, dig=False, dict_rep=False, certify=False, |
1356 | | verbosity=0, use_indicator_function=True, sparse=False, |
1357 | | base=False, order=False): |
1358 | | """ |
1359 | | Assumes that the vertex set of G is 0,1,...,n-1 for some n. |
1360 | | |
1361 | | Note that this conflicts with the SymmetricGroup we are using to |
1362 | | represent automorphisms. The solution is to let the group act on |
1363 | | the set 1,2,...,n, under the assumption n = 0. |
1364 | | |
1365 | | INPUT: lab- if True, return the canonical label in addition to the |
1366 | | auto- morphism group. dig- if True, does not use Lemma 2.25 in [1], |
1367 | | and the algorithm is valid for digraphs and graphs with loops. |
1368 | | dict_rep- if True, explain which vertices are which elements of |
1369 | | the set 1,2,...,n in the representation of the automorphism group. |
1370 | | certify- if True, return the relabeling from G to its canonical |
1371 | | label. Forces lab=True. verbosity- 0 - print nothing 1 - display |
1372 | | state trace 2 - with timings 3 - display partition nests 4 - |
1373 | | display orbit partition 5 - plot the part of the tree traversed |
1374 | | during search |
1375 | | |
1376 | | |
1377 | | - ``use_indicator_function`` - option to turn off |
1378 | | indicator function (False - slower) |
1379 | | |
1380 | | - ``sparse`` - whether to use sparse or dense |
1381 | | representation of the graph (ignored if G is already a CGraph - see |
1382 | | sage.graphs.base) |
1383 | | |
1384 | | - ``base`` - whether to return the first sequence of |
1385 | | split vertices (used in computing the order of the group) |
1386 | | |
1387 | | - ``order`` - whether to return the order of the |
1388 | | automorphism group |
1389 | | |
1390 | | |
1391 | | STATE DIAGRAM:: |
1392 | | |
1393 | | sage: SD = DiGraph( { 1:[18,2], 2:[5,3], 3:[4,6], 4:[7,2], 5:[4], 6:[13,12], 7:[18,8,10], 8:[6,9,10], 9:[6], 10:[11,13], 11:[12], 12:[13], 13:[17,14], 14:[16,15], 15:[2], 16:[13], 17:[15,13], 18:[13] }, implementation='networkx' ) |
1394 | | sage: SD.set_edge_label(1, 18, 'discrete') |
1395 | | sage: SD.set_edge_label(4, 7, 'discrete') |
1396 | | sage: SD.set_edge_label(2, 5, 'h = 0') |
1397 | | sage: SD.set_edge_label(7, 18, 'h = 0') |
1398 | | sage: SD.set_edge_label(7, 10, 'aut') |
1399 | | sage: SD.set_edge_label(8, 10, 'aut') |
1400 | | sage: SD.set_edge_label(8, 9, 'label') |
1401 | | sage: SD.set_edge_label(8, 6, 'no label') |
1402 | | sage: SD.set_edge_label(13, 17, 'k > h') |
1403 | | sage: SD.set_edge_label(13, 14, 'k = h') |
1404 | | sage: SD.set_edge_label(17, 15, 'v_k finite') |
1405 | | sage: SD.set_edge_label(14, 15, 'v_k m.c.r.') |
1406 | | sage: posn = {1:[ 3,-3], 2:[0,2], 3:[0, 13], 4:[3,9], 5:[3,3], 6:[16, 13], 7:[6,1], 8:[6,6], 9:[6,11], 10:[9,1], 11:[10,6], 12:[13,6], 13:[16,2], 14:[10,-6], 15:[0,-10], 16:[14,-6], 17:[16,-10], 18:[6,-4]} |
1407 | | sage: SD.plot(pos=posn, vertex_size=400, vertex_colors={'#FFFFFF':range(1,19)}, edge_labels=True) |
1408 | | |
1409 | | .. note:: |
1410 | | |
1411 | | There is a function, called test_refine, that has the same |
1412 | | signature as _refine. It calls _refine, then checks to make |
1413 | | sure the output is sane. To use this, simply add 'test' to the |
1414 | | two places this algorithm calls the function (states 1 and 2). |
1415 | | |
1416 | | EXAMPLES: The following example is due to Chris Godsil:: |
1417 | | |
1418 | | sage: HS = graphs.HoffmanSingletonGraph() |
1419 | | sage: clqs = (HS.complement()).cliques() |
1420 | | sage: alqs = [Set(c) for c in clqs if len(c) == 15] |
1421 | | sage: Y = Graph([alqs, lambda s,t: len(s.intersection(t))==0], implementation='networkx') |
1422 | | sage: Y0,Y1 = Y.connected_components_subgraphs() |
1423 | | sage: Y0.is_isomorphic(Y1) |
1424 | | True |
1425 | | sage: Y0.is_isomorphic(HS) |
1426 | | True |
1427 | | |
1428 | | :: |
1429 | | |
1430 | | sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree |
1431 | | sage: from sage.graphs.base.sparse_graph import SparseGraph |
1432 | | sage: from sage.graphs.base.dense_graph import DenseGraph |
1433 | | sage: from sage.groups.perm_gps.permgroup import PermutationGroup |
1434 | | sage: from sage.graphs.graph_isom import perm_group_elt |
1435 | | |
1436 | | :: |
1437 | | |
1438 | | sage: G = graphs.DodecahedralGraph() |
1439 | | sage: GD = DenseGraph(20) |
1440 | | sage: GS = SparseGraph(20) |
1441 | | sage: for i,j,_ in G.edge_iterator(): |
1442 | | ... GD.add_arc(i,j); GD.add_arc(j,i) |
1443 | | ... GS.add_arc(i,j); GS.add_arc(j,i) |
1444 | | sage: Pi=[range(20)] |
1445 | | sage: a,b = search_tree(G, Pi) |
1446 | | sage: asp,bsp = search_tree(GS, Pi) |
1447 | | sage: ade,bde = search_tree(GD, Pi) |
1448 | | sage: bsg = Graph(implementation='networkx') |
1449 | | sage: bdg = Graph(implementation='networkx') |
1450 | | sage: for i in range(20): |
1451 | | ... for j in range(20): |
1452 | | ... if bsp.has_arc(i,j): |
1453 | | ... bsg.add_edge(i,j) |
1454 | | ... if bde.has_arc(i,j): |
1455 | | ... bdg.add_edge(i,j) |
1456 | | sage: print a, b.graph6_string() |
1457 | | [[0, 19, 3, 2, 6, 5, 4, 17, 18, 11, 10, 9, 13, 12, 16, 15, 14, 7, 8, 1], [0, 1, 8, 9, 13, 14, 7, 6, 2, 3, 19, 18, 17, 4, 5, 15, 16, 12, 11, 10], [1, 8, 9, 10, 11, 12, 13, 14, 7, 6, 2, 3, 4, 5, 15, 16, 17, 18, 19, 0]] S?[PG__OQ@?_?_?P?CO?_?AE?EC?Ac?@O |
1458 | | sage: a == asp |
1459 | | True |
1460 | | sage: a == ade |
1461 | | True |
1462 | | sage: b == bsg |
1463 | | True |
1464 | | sage: b == bdg |
1465 | | True |
1466 | | sage: c = search_tree(G, Pi, lab=False) |
1467 | | sage: print c |
1468 | | [[0, 19, 3, 2, 6, 5, 4, 17, 18, 11, 10, 9, 13, 12, 16, 15, 14, 7, 8, 1], [0, 1, 8, 9, 13, 14, 7, 6, 2, 3, 19, 18, 17, 4, 5, 15, 16, 12, 11, 10], [1, 8, 9, 10, 11, 12, 13, 14, 7, 6, 2, 3, 4, 5, 15, 16, 17, 18, 19, 0]] |
1469 | | sage: DodecAut = PermutationGroup([perm_group_elt(aa) for aa in a]) |
1470 | | sage: DodecAut.character_table() |
1471 | | [ 1 1 1 1 1 1 1 1 1 1] |
1472 | | [ 1 -1 1 1 -1 1 -1 1 -1 -1] |
1473 | | [ 3 -1 0 -1 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 3] |
1474 | | [ 3 -1 0 -1 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 3] |
1475 | | [ 3 1 0 -1 zeta5^3 + zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 -zeta5^3 - zeta5^2 - 1 -3] |
1476 | | [ 3 1 0 -1 -zeta5^3 - zeta5^2 - 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 zeta5^3 + zeta5^2 -3] |
1477 | | [ 4 0 1 0 -1 -1 1 -1 -1 4] |
1478 | | [ 4 0 1 0 1 -1 -1 -1 1 -4] |
1479 | | [ 5 1 -1 1 0 0 -1 0 0 5] |
1480 | | [ 5 -1 -1 1 0 0 1 0 0 -5] |
1481 | | |
1482 | | sage: DodecAut2 = PermutationGroup([perm_group_elt(cc) for cc in c]) |
1483 | | sage: DodecAut2.character_table() |
1484 | | [ 1 1 1 1 1 1 1 1 1 1] |
1485 | | [ 1 -1 1 1 -1 1 -1 1 -1 -1] |
1486 | | [ 3 -1 0 -1 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 3] |
1487 | | [ 3 -1 0 -1 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 3] |
1488 | | [ 3 1 0 -1 zeta5^3 + zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 -zeta5^3 - zeta5^2 - 1 -3] |
1489 | | [ 3 1 0 -1 -zeta5^3 - zeta5^2 - 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 zeta5^3 + zeta5^2 -3] |
1490 | | [ 4 0 1 0 -1 -1 1 -1 -1 4] |
1491 | | [ 4 0 1 0 1 -1 -1 -1 1 -4] |
1492 | | [ 5 1 -1 1 0 0 -1 0 0 5] |
1493 | | [ 5 -1 -1 1 0 0 1 0 0 -5] |
1494 | | |
1495 | | |
1496 | | :: |
1497 | | |
1498 | | sage: G = graphs.PetersenGraph() |
1499 | | sage: Pi=[range(10)] |
1500 | | sage: a,b = search_tree(G, Pi) |
1501 | | sage: print a, b.graph6_string() |
1502 | | [[0, 1, 2, 7, 5, 4, 6, 3, 9, 8], [0, 1, 6, 8, 5, 4, 2, 9, 3, 7], [0, 4, 3, 8, 5, 1, 9, 2, 6, 7], [1, 0, 4, 9, 6, 2, 5, 3, 7, 8]] I@OZCMgs? |
1503 | | sage: c = search_tree(G, Pi, lab=False) |
1504 | | sage: PAut = PermutationGroup([perm_group_elt(aa) for aa in a]) |
1505 | | sage: PAut.character_table() |
1506 | | [ 1 1 1 1 1 1 1] |
1507 | | [ 1 -1 1 -1 1 -1 1] |
1508 | | [ 4 -2 0 1 1 0 -1] |
1509 | | [ 4 2 0 -1 1 0 -1] |
1510 | | [ 5 1 1 1 -1 -1 0] |
1511 | | [ 5 -1 1 -1 -1 1 0] |
1512 | | [ 6 0 -2 0 0 0 1] |
1513 | | sage: PAut = PermutationGroup([perm_group_elt(cc) for cc in c]) |
1514 | | sage: PAut.character_table() |
1515 | | [ 1 1 1 1 1 1 1] |
1516 | | [ 1 -1 1 -1 1 -1 1] |
1517 | | [ 4 -2 0 1 1 0 -1] |
1518 | | [ 4 2 0 -1 1 0 -1] |
1519 | | [ 5 1 1 1 -1 -1 0] |
1520 | | [ 5 -1 1 -1 -1 1 0] |
1521 | | [ 6 0 -2 0 0 0 1] |
1522 | | |
1523 | | :: |
1524 | | |
1525 | | sage: G = graphs.CubeGraph(3) |
1526 | | sage: Pi = [] |
1527 | | sage: for i in range(8): |
1528 | | ... b = Integer(i).binary() |
1529 | | ... Pi.append(b.zfill(3)) |
1530 | | ... |
1531 | | sage: Pi = [Pi] |
1532 | | sage: a,b = search_tree(G, Pi) |
1533 | | sage: print a, b.graph6_string() |
1534 | | [[0, 2, 1, 3, 4, 6, 5, 7], [0, 1, 4, 5, 2, 3, 6, 7], [1, 0, 3, 2, 5, 4, 7, 6]] GIQ\T_ |
1535 | | sage: c = search_tree(G, Pi, lab=False) |
1536 | | |
1537 | | :: |
1538 | | |
1539 | | sage: PermutationGroup([perm_group_elt(aa) for aa in a]).order() |
1540 | | 48 |
1541 | | sage: PermutationGroup([perm_group_elt(cc) for cc in c]).order() |
1542 | | 48 |
1543 | | sage: DodecAut.order() |
1544 | | 120 |
1545 | | sage: PAut.order() |
1546 | | 120 |
1547 | | |
1548 | | :: |
1549 | | |
1550 | | sage: D = graphs.DodecahedralGraph() |
1551 | | sage: a,b,c = search_tree(D, [range(20)], certify=True) |
1552 | | sage: from sage.plot.plot import GraphicsArray |
1553 | | sage: from sage.graphs.graph_fast import spring_layout_fast |
1554 | | sage: position_D = spring_layout_fast(D) |
1555 | | sage: position_b = {} |
1556 | | sage: for vert in position_D: |
1557 | | ... position_b[c[vert]] = position_D[vert] |
1558 | | sage: graphics_array([D.plot(pos=position_D), b.plot(pos=position_b)]).show() |
1559 | | sage: c |
1560 | | {0: 0, 1: 19, 2: 16, 3: 15, 4: 9, 5: 1, 6: 10, 7: 8, 8: 14, 9: 12, 10: 17, 11: 11, 12: 5, 13: 6, 14: 2, 15: 4, 16: 3, 17: 7, 18: 13, 19: 18} |
1561 | | |
1562 | | BENCHMARKS: The following examples are given to check modifications |
1563 | | to the algorithm for optimization. |
1564 | | |
1565 | | :: |
1566 | | |
1567 | | sage: G = Graph({0:[]}) |
1568 | | sage: Pi = [[0]] |
1569 | | sage: a,b = search_tree(G, Pi) |
1570 | | sage: print a, b.graph6_string() |
1571 | | [] @ |
1572 | | sage: a,b = search_tree(G, Pi, dig=True) |
1573 | | sage: print a, b.graph6_string() |
1574 | | [] @ |
1575 | | sage: search_tree(G, Pi, lab=False) |
1576 | | [] |
1577 | | |
1578 | | :: |
1579 | | |
1580 | | sage: from sage.graphs.graph_isom import all_labeled_graphs, all_ordered_partitions |
1581 | | |
1582 | | :: |
1583 | | |
1584 | | sage: graph2 = all_labeled_graphs(2) |
1585 | | sage: part2 = all_ordered_partitions(range(2)) |
1586 | | sage: for G in graph2: |
1587 | | ... for Pi in part2: |
1588 | | ... a,b = search_tree(G, Pi) |
1589 | | ... c,d = search_tree(G, Pi, dig=True) |
1590 | | ... e = search_tree(G, Pi, lab=False) |
1591 | | ... a = str(a); b = b.graph6_string(); c = str(c); d = d.graph6_string(); e = str(e) |
1592 | | ... print a.ljust(15), b.ljust(5), c.ljust(15), d.ljust(5), e.ljust(15) |
1593 | | [] A? [] A? [] |
1594 | | [] A? [] A? [] |
1595 | | [[1, 0]] A? [[1, 0]] A? [[1, 0]] |
1596 | | [[1, 0]] A? [[1, 0]] A? [[1, 0]] |
1597 | | [] A_ [] A_ [] |
1598 | | [] A_ [] A_ [] |
1599 | | [[1, 0]] A_ [[1, 0]] A_ [[1, 0]] |
1600 | | [[1, 0]] A_ [[1, 0]] A_ [[1, 0]] |
1601 | | |
1602 | | :: |
1603 | | |
1604 | | sage: graph3 = all_labeled_graphs(3) |
1605 | | sage: part3 = all_ordered_partitions(range(3)) |
1606 | | sage: for G in graph3: |
1607 | | ... for Pi in part3: |
1608 | | ... a,b = search_tree(G, Pi) |
1609 | | ... c,d = search_tree(G, Pi, dig=True) |
1610 | | ... e = search_tree(G, Pi, lab=False) |
1611 | | ... a = str(a); b = b.graph6_string(); c = str(c); d = d.graph6_string(); e = str(e) |
1612 | | ... print a.ljust(15), b.ljust(5), c.ljust(15), d.ljust(5), e.ljust(15) |
1613 | | [] B? [] B? [] |
1614 | | [] B? [] B? [] |
1615 | | [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]] |
1616 | | [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]] |
1617 | | [] B? [] B? [] |
1618 | | [] B? [] B? [] |
1619 | | [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]] |
1620 | | [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]] |
1621 | | [] B? [] B? [] |
1622 | | [] B? [] B? [] |
1623 | | [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]] |
1624 | | [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]] |
1625 | | [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]] |
1626 | | [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]] |
1627 | | [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]] |
1628 | | [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]] |
1629 | | [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]] |
1630 | | [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]] |
1631 | | [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] |
1632 | | [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] |
1633 | | [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] |
1634 | | [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] |
1635 | | [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] |
1636 | | [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] |
1637 | | [] BG [] BG [] |
1638 | | [] BG [] BG [] |
1639 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1640 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1641 | | [] BO [] BO [] |
1642 | | [] B_ [] B_ [] |
1643 | | [] BO [] BO [] |
1644 | | [] BO [] BO [] |
1645 | | [] BO [] BO [] |
1646 | | [] B_ [] B_ [] |
1647 | | [] BO [] BO [] |
1648 | | [] BO [] BO [] |
1649 | | [] BG [] BG [] |
1650 | | [] BG [] BG [] |
1651 | | [] BG [] BG [] |
1652 | | [[0, 2, 1]] B_ [[0, 2, 1]] B_ [[0, 2, 1]] |
1653 | | [] BG [] BG [] |
1654 | | [[0, 2, 1]] B_ [[0, 2, 1]] B_ [[0, 2, 1]] |
1655 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1656 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1657 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1658 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1659 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1660 | | [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]] |
1661 | | [] BO [] BO [] |
1662 | | [] B_ [] B_ [] |
1663 | | [] BO [] BO [] |
1664 | | [] BO [] BO [] |
1665 | | [] BG [] BG [] |
1666 | | [] BG [] BG [] |
1667 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1668 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1669 | | [] B_ [] B_ [] |
1670 | | [] BO [] BO [] |
1671 | | [] BO [] BO [] |
1672 | | [] BO [] BO [] |
1673 | | [] BG [] BG [] |
1674 | | [[2, 1, 0]] B_ [[2, 1, 0]] B_ [[2, 1, 0]] |
1675 | | [] BG [] BG [] |
1676 | | [] BG [] BG [] |
1677 | | [[2, 1, 0]] B_ [[2, 1, 0]] B_ [[2, 1, 0]] |
1678 | | [] BG [] BG [] |
1679 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1680 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1681 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1682 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1683 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1684 | | [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]] |
1685 | | [] BW [] BW [] |
1686 | | [] Bg [] Bg [] |
1687 | | [] BW [] BW [] |
1688 | | [] BW [] BW [] |
1689 | | [] BW [] BW [] |
1690 | | [] Bg [] Bg [] |
1691 | | [] BW [] BW [] |
1692 | | [] BW [] BW [] |
1693 | | [] Bo [] Bo [] |
1694 | | [] Bo [] Bo [] |
1695 | | [[1, 0, 2]] Bo [[1, 0, 2]] Bo [[1, 0, 2]] |
1696 | | [[1, 0, 2]] Bo [[1, 0, 2]] Bo [[1, 0, 2]] |
1697 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1698 | | [] Bg [] Bg [] |
1699 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1700 | | [] Bg [] Bg [] |
1701 | | [] Bg [] Bg [] |
1702 | | [] Bg [] Bg [] |
1703 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1704 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1705 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1706 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1707 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1708 | | [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]] |
1709 | | [] B_ [] B_ [] |
1710 | | [] BO [] BO [] |
1711 | | [] BO [] BO [] |
1712 | | [] BO [] BO [] |
1713 | | [] B_ [] B_ [] |
1714 | | [] BO [] BO [] |
1715 | | [] BO [] BO [] |
1716 | | [] BO [] BO [] |
1717 | | [] BG [] BG [] |
1718 | | [] BG [] BG [] |
1719 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1720 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1721 | | [[1, 0, 2]] B_ [[1, 0, 2]] B_ [[1, 0, 2]] |
1722 | | [] BG [] BG [] |
1723 | | [[1, 0, 2]] B_ [[1, 0, 2]] B_ [[1, 0, 2]] |
1724 | | [] BG [] BG [] |
1725 | | [] BG [] BG [] |
1726 | | [] BG [] BG [] |
1727 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1728 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1729 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1730 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1731 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1732 | | [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]] |
1733 | | [] Bg [] Bg [] |
1734 | | [] BW [] BW [] |
1735 | | [] BW [] BW [] |
1736 | | [] BW [] BW [] |
1737 | | [] Bo [] Bo [] |
1738 | | [] Bo [] Bo [] |
1739 | | [[2, 1, 0]] Bo [[2, 1, 0]] Bo [[2, 1, 0]] |
1740 | | [[2, 1, 0]] Bo [[2, 1, 0]] Bo [[2, 1, 0]] |
1741 | | [] BW [] BW [] |
1742 | | [] Bg [] Bg [] |
1743 | | [] BW [] BW [] |
1744 | | [] BW [] BW [] |
1745 | | [] Bg [] Bg [] |
1746 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1747 | | [] Bg [] Bg [] |
1748 | | [] Bg [] Bg [] |
1749 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1750 | | [] Bg [] Bg [] |
1751 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1752 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1753 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1754 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1755 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1756 | | [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]] |
1757 | | [] Bo [] Bo [] |
1758 | | [] Bo [] Bo [] |
1759 | | [[0, 2, 1]] Bo [[0, 2, 1]] Bo [[0, 2, 1]] |
1760 | | [[0, 2, 1]] Bo [[0, 2, 1]] Bo [[0, 2, 1]] |
1761 | | [] Bg [] Bg [] |
1762 | | [] BW [] BW [] |
1763 | | [] BW [] BW [] |
1764 | | [] BW [] BW [] |
1765 | | [] Bg [] Bg [] |
1766 | | [] BW [] BW [] |
1767 | | [] BW [] BW [] |
1768 | | [] BW [] BW [] |
1769 | | [] Bg [] Bg [] |
1770 | | [] Bg [] Bg [] |
1771 | | [] Bg [] Bg [] |
1772 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1773 | | [] Bg [] Bg [] |
1774 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1775 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1776 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1777 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1778 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1779 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1780 | | [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]] |
1781 | | [] Bw [] Bw [] |
1782 | | [] Bw [] Bw [] |
1783 | | [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]] |
1784 | | [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]] |
1785 | | [] Bw [] Bw [] |
1786 | | [] Bw [] Bw [] |
1787 | | [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]] |
1788 | | [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]] |
1789 | | [] Bw [] Bw [] |
1790 | | [] Bw [] Bw [] |
1791 | | [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]] |
1792 | | [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]] |
1793 | | [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]] |
1794 | | [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]] |
1795 | | [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]] |
1796 | | [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]] |
1797 | | [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]] |
1798 | | [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]] |
1799 | | [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] |
1800 | | [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] |
1801 | | [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] |
1802 | | [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] |
1803 | | [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] |
1804 | | [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] |
1805 | | |
1806 | | :: |
1807 | | |
1808 | | sage: C = graphs.CubeGraph(1) |
1809 | | sage: gens = search_tree(C, [C.vertices()], lab=False) |
1810 | | sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() |
1811 | | 2 |
1812 | | sage: C = graphs.CubeGraph(2) |
1813 | | sage: gens = search_tree(C, [C.vertices()], lab=False) |
1814 | | sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() |
1815 | | 8 |
1816 | | sage: C = graphs.CubeGraph(3) |
1817 | | sage: gens = search_tree(C, [C.vertices()], lab=False) |
1818 | | sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() |
1819 | | 48 |
1820 | | sage: C = graphs.CubeGraph(4) |
1821 | | sage: gens = search_tree(C, [C.vertices()], lab=False) |
1822 | | sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() |
1823 | | 384 |
1824 | | sage: C = graphs.CubeGraph(5) |
1825 | | sage: gens = search_tree(C, [C.vertices()], lab=False) |
1826 | | sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() |
1827 | | 3840 |
1828 | | sage: C = graphs.CubeGraph(6) |
1829 | | sage: gens = search_tree(C, [C.vertices()], lab=False) |
1830 | | sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() |
1831 | | 46080 |
1832 | | |
1833 | | One can also turn off the indicator function (note- this will take |
1834 | | longer) |
1835 | | |
1836 | | :: |
1837 | | |
1838 | | sage: D1 = DiGraph({0:[2],2:[0],1:[1]}, loops=True) |
1839 | | sage: D2 = DiGraph({1:[2],2:[1],0:[0]}, loops=True) |
1840 | | sage: a,b = search_tree(D1, [D1.vertices()], use_indicator_function=False) |
1841 | | sage: c,d = search_tree(D2, [D2.vertices()], use_indicator_function=False) |
1842 | | sage: b==d |
1843 | | True |
1844 | | |
1845 | | Previously a bug, now the output is correct:: |
1846 | | |
1847 | | sage: G = Graph('^????????????????????{??N??@w??FaGa?PCO@CP?AGa?_QO?Q@G?CcA??cc????Bo????{????F_') |
1848 | | sage: perm = {3:15, 15:3} |
1849 | | sage: H = G.relabel(perm, inplace=False) |
1850 | | sage: G.canonical_label() == H.canonical_label() |
1851 | | True |
1852 | | |
1853 | | Another former bug:: |
1854 | | |
1855 | | sage: Graph('Fll^G').canonical_label() |
1856 | | Graph on 7 vertices |
1857 | | |
1858 | | :: |
1859 | | |
1860 | | sage: g = Graph(21) |
1861 | | sage: g.automorphism_group(return_group=False, order=True) |
1862 | | 51090942171709440000 |
1863 | | """ |
1864 | | cdef int i, j, m # local variables |
1865 | | cdef int uif = 1 if use_indicator_function else 0 |
1866 | | cdef int _base = 1 if base else 0 |
1867 | | |
1868 | | cdef OrbitPartition Theta |
1869 | | cdef int index = 0 # see Theorem 2.33 in [1] |
1870 | | size = Integer(1) |
1871 | | |
1872 | | cdef int L = 100 # memory limit for storing values from fix and mcr: |
1873 | | # Phi and Omega store specific information about some |
1874 | | # of the automorphisms we already know about, and they |
1875 | | # are arrays of length L |
1876 | | cdef int **Phi # stores the fixed point sets of each automorphism |
1877 | | cdef int **Omega # stores the minimal elements of each cell of the |
1878 | | # orbit partition |
1879 | | cdef int l = -1 # current index for storing values in Phi and Omega- |
1880 | | # we start at -1 so that when we increment first, |
1881 | | # the first place we write to is 0. |
1882 | | cdef int **W # for each k, W[k] is a list of the vertices to be searched |
1883 | | # down from the current partition nest, at k |
1884 | | # Phi and Omega are ultimately used to make the size of W |
1885 | | # as small as possible |
1886 | | |
1887 | | cdef PartitionStack nu, zeta, rho |
1888 | | cdef int k_rho # the number of partitions in rho |
1889 | | cdef int h = -1 # longest common ancestor of zeta and nu: |
1890 | | # zeta[h] == nu[h], zeta[h+1] != nu[h+1] |
1891 | | cdef int hb # longest common ancestor of rho and nu: |
1892 | | # rho[hb] == nu[hb], rho[hb+1] != nu[hb+1] |
1893 | | cdef int hh = 1 # the height of the oldest ancestor of nu |
1894 | | # satisfying Lemma 2.25 in [1] |
1895 | | cdef int ht # smallest such that all descendants of zeta[ht] |
1896 | | # are known to be equivalent |
1897 | | |
1898 | | cdef mpz_t *Lambda_mpz, *zf_mpz, *zb_mpz # for tracking indicator values |
1899 | | # zf and zb are indicator vectors remembering Lambda[k] for zeta and rho, |
1900 | | # respectively |
1901 | | cdef int hzf # the max height for which Lambda and zf agree |
1902 | | cdef int hzb = -1 # the max height for which Lambda and zb agree |
1903 | | |
1904 | | cdef CGraph M |
1905 | | cdef int *gamma # for storing permutations |
1906 | | cdef int *alpha # for storing pointers to cells of nu[k]: |
1907 | | # allocated to be length 4*n + 1 for scratch (see functions |
1908 | | # _sort_by_function and _refine) |
1909 | | cdef int *v # list of vertices determining nu |
1910 | | cdef int *e # 0 or 1, see states 12 and 17 |
1911 | | cdef int state # keeps track of place in algorithm |
1912 | | cdef int _dig, tvh, n |
1913 | | |
1914 | | if isinstance(G, GenericGraph): |
1915 | | n = G.order() |
1916 | | elif isinstance(G, CGraph): |
1917 | | M = G |
1918 | | n = M.num_verts |
1919 | | else: |
1920 | | raise TypeError("G must be a Sage graph.") |
1921 | | |
1922 | | # trivial case |
1923 | | if n == 0: |
1924 | | if not (lab or dict_rep or certify or base or order): |
1925 | | return [[]] |
1926 | | output_tuple = [[[]]] |
1927 | | if dict_rep: |
1928 | | output_tuple.append({}) |
1929 | | if lab: |
1930 | | output_tuple.append(G.copy()) |
1931 | | if certify: |
1932 | | output_tuple.append({}) |
1933 | | if base: |
1934 | | output_tuple.append([]) |
1935 | | if order: |
1936 | | output_tuple.append(1) |
1937 | | return tuple(output_tuple) |
1938 | | |
1939 | | # allocate int pointers |
1940 | | W = <int **> sage_malloc( n * sizeof(int *) ) |
1941 | | Phi = <int **> sage_malloc( L * sizeof(int *) ) |
1942 | | Omega = <int **> sage_malloc( L * sizeof(int *) ) |
1943 | | |
1944 | | # allocate GMP int pointers |
1945 | | Lambda_mpz = <mpz_t *> sage_malloc( (n+2) * sizeof(mpz_t) ) |
1946 | | zf_mpz = <mpz_t *> sage_malloc( (n+2) * sizeof(mpz_t) ) |
1947 | | zb_mpz = <mpz_t *> sage_malloc( (n+2) * sizeof(mpz_t) ) |
1948 | | |
1949 | | # check for memory errors |
1950 | | if not (W and Phi and Omega and Lambda_mpz and zf_mpz and zb_mpz): |
1951 | | sage_free(Lambda_mpz) |
1952 | | sage_free(zf_mpz) |
1953 | | sage_free(zb_mpz) |
1954 | | sage_free(W) |
1955 | | sage_free(Phi) |
1956 | | sage_free(Omega) |
1957 | | raise MemoryError("Error allocating memory.") |
1958 | | |
1959 | | # allocate int arrays |
1960 | | gamma = <int *> sage_malloc( n * sizeof(int) ) |
1961 | | W[0] = <int *> sage_malloc( (n*n) * sizeof(int) ) |
1962 | | Phi[0] = <int *> sage_malloc( (L*n) * sizeof(int) ) |
1963 | | Omega[0] = <int *> sage_malloc( (L*n) * sizeof(int) ) |
1964 | | alpha = <int *> sage_malloc( (4*n + 1) * sizeof(int) ) |
1965 | | v = <int *> sage_malloc( n * sizeof(int) ) |
1966 | | e = <int *> sage_malloc( n * sizeof(int) ) |
1967 | | |
1968 | | # check for memory errors |
1969 | | if not (gamma and W[0] and Phi[0] and Omega[0] and alpha and v and e): |
1970 | | sage_free(gamma) |
1971 | | sage_free(W[0]) |
1972 | | sage_free(Phi[0]) |
1973 | | sage_free(Omega[0]) |
1974 | | sage_free(alpha) |
1975 | | sage_free(v) |
1976 | | sage_free(e) |
1977 | | sage_free(Lambda_mpz) |
1978 | | sage_free(zf_mpz) |
1979 | | sage_free(zb_mpz) |
1980 | | sage_free(W) |
1981 | | sage_free(Phi) |
1982 | | sage_free(Omega) |
1983 | | raise MemoryError("Error allocating memory.") |
1984 | | |
1985 | | # setup double index arrays |
1986 | | for i from 0 < i < n: |
1987 | | W[i] = W[0] + n*i |
1988 | | for i from 0 < i < L: |
1989 | | Phi[i] = Phi[0] + n*i |
1990 | | for i from 0 < i < L: |
1991 | | Omega[i] = Omega[0] + n*i |
1992 | | |
1993 | | # allocate GMP ints |
1994 | | for i from 0 <= i < n+2: |
1995 | | mpz_init(Lambda_mpz[i]) |
1996 | | mpz_init_set_si(zf_mpz[i], -1) # correspond to default values of |
1997 | | mpz_init_set_si(zb_mpz[i], -1) # "infinity" |
1998 | | # Note that there is a potential memory leak here - if a particular |
1999 | | # mpz fails to allocate, this is not checked for |
2000 | | |
2001 | | if isinstance(G, GenericGraph): |
2002 | | # relabel vertices to the set {0,...,n-1} |
2003 | | G = G.copy() |
2004 | | ffrom = G.relabel(return_map=True) |
2005 | | to = {} |
2006 | | for vvv in ffrom.iterkeys(): |
2007 | | to[ffrom[vvv]] = vvv |
2008 | | Pi2 = [] |
2009 | | for cell in Pi: |
2010 | | Pi2.append([ffrom[c] for c in cell]) |
2011 | | Pi = Pi2 |
2012 | | if sparse: |
2013 | | M = SparseGraph(n) |
2014 | | else: |
2015 | | M = DenseGraph(n) |
2016 | | if isinstance(G, Graph): |
2017 | | for i, j, la in G.edge_iterator(): |
2018 | | M.add_arc_unsafe(i,j) |
2019 | | M.add_arc_unsafe(j,i) |
2020 | | elif isinstance(G, DiGraph): |
2021 | | for i, j, la in G.edge_iterator(): |
2022 | | M.add_arc_unsafe(i,j) |
2023 | | |
2024 | | # initialize W |
2025 | | for i from 0 <= i < n: |
2026 | | for j from 0 <= j < n: |
2027 | | W[i][j] = 0 |
2028 | | |
2029 | | # set up the rest of the variables |
2030 | | nu = PartitionStack(Pi) |
2031 | | Theta = OrbitPartition(n) |
2032 | | output = [] |
2033 | | if dig: _dig = 1 |
2034 | | else: _dig = 0 |
2035 | | if certify: |
2036 | | lab=True |
2037 | | if _base: |
2038 | | base_set = [] |
2039 | | |
2040 | | if verbosity > 1: |
2041 | | t = cputime() |
2042 | | if verbosity > 2: |
2043 | | rho = PartitionStack(n) |
2044 | | zeta = PartitionStack(n) |
2045 | | state = 1 |
2046 | | while state != -1: |
2047 | | if verbosity > 0: |
2048 | | print '-----' |
2049 | | print 'k: ' + str(nu.k) |
2050 | | print 'k_rho: ' + str(k_rho) |
2051 | | print 'hh', hh |
2052 | | print 'nu:' |
2053 | | print nu |
2054 | | if verbosity >= 2: |
2055 | | t = cputime(t) |
2056 | | print 'time:', t |
2057 | | if verbosity >= 3: |
2058 | | print 'zeta:' |
2059 | | print [zeta.entries[iii] for iii in range(n)] |
2060 | | print [zeta.levels[iii] for iii in range(n)] |
2061 | | print 'rho' |
2062 | | print [rho.entries[iii] for iii in range(n)] |
2063 | | print [rho.levels[iii] for iii in range(n)] |
2064 | | if verbosity >= 4: |
2065 | | Thetarep = [] |
2066 | | for i from 0 <= i < n: |
2067 | | j = Theta._find(i) |
2068 | | didit = False |
2069 | | for celll in Thetarep: |
2070 | | if celll[0] == j: |
2071 | | celll.append(i) |
2072 | | didit = True |
2073 | | if not didit: |
2074 | | Thetarep.append([j]) |
2075 | | print 'Theta: ', str(Thetarep) |
2076 | | if verbosity >= 5: |
2077 | | if state == 1: |
2078 | | verbose_first_time = True |
2079 | | verbose_just_refined = False |
2080 | | elif verbose_first_time: |
2081 | | verbose_first_time = False |
2082 | | # here we have just gone through step 1, and must now begin |
2083 | | # to record information about the tree |
2084 | | ST_vis = DiGraph() |
2085 | | ST_vis_heights = {0:[nu.repr_at_k(0)]} |
2086 | | ST_vis.add_vertex(nu.repr_at_k(0)) |
2087 | | ST_vis_current_vertex = nu.repr_at_k(0) |
2088 | | ST_vis_current_level = 0 |
2089 | | #ST_vis.show(vertex_size=0) |
2090 | | if state == 2: |
2091 | | verbose_just_refined = True |
2092 | | elif verbose_just_refined: |
2093 | | verbose_just_refined = False |
2094 | | # here we have gone through step 2, and must record the |
2095 | | # refinement just made |
2096 | | while ST_vis_current_level > nu.k-1: |
2097 | | ST_vis_current_vertex = ST_vis.predecessors(ST_vis_current_vertex)[0] |
2098 | | ST_vis_current_level -= 1 |
2099 | | if ST_vis_heights.has_key(nu.k): |
2100 | | ST_vis_heights[nu.k].append(nu.repr_at_k(nu.k)) |
2101 | | else: |
2102 | | ST_vis_heights[nu.k] = [nu.repr_at_k(nu.k)] |
2103 | | ST_vis.add_edge(ST_vis_current_vertex, nu.repr_at_k(nu.k), '%d,%d'%(ST_vis_splitvert, ST_vis_invariant)) |
2104 | | ST_vis_current_vertex = nu.repr_at_k(nu.k) |
2105 | | ST_vis_current_level += 1 |
2106 | | if state == 13 and nu.k == -1: |
2107 | | ST_vis_new_heights = {} |
2108 | | for ST_vis_k in ST_vis_heights: |
2109 | | ST_vis_new_heights[-ST_vis_k] = ST_vis_heights[ST_vis_k] |
2110 | | ST_vis.show(vertex_size=0, heights=ST_vis_new_heights, figsize=[30,10], edge_labels=True, edge_colors={(.6,.6,.6):ST_vis.edges()}) |
2111 | | print '-----' |
2112 | | print 'state:', state |
2113 | | |
2114 | | |
2115 | | if state == 1: # Entry point to algorithm |
2116 | | # get alpha to point to cells of nu |
2117 | | j = 1 |
2118 | | alpha[0] = 0 |
2119 | | for i from 0 < i < n: |
2120 | | if nu.levels[i-1] == 0: |
2121 | | alpha[j] = i |
2122 | | j += 1 |
2123 | | alpha[j] = -1 |
2124 | | |
2125 | | # "nu[0] := R(G, Pi, Pi)" |
2126 | | nu._refine(alpha, n, M, _dig, uif) |
2127 | | |
2128 | | if not _dig: |
2129 | | if nu._sat_225(n): hh = nu.k |
2130 | | if nu._is_discrete(): state = 18; continue |
2131 | | |
2132 | | # store the first smallest nontrivial cell in W[k], and set v[k] |
2133 | | # equal to its minimum element |
2134 | | v[nu.k] = nu._first_smallest_nontrivial(W[nu.k], n) |
2135 | | mpz_set_ui(Lambda_mpz[nu.k], 0) |
2136 | | e[nu.k] = 0 # see state 12, and 17 |
2137 | | state = 2 |
2138 | | |
2139 | | elif state == 2: # Move down the search tree one level by refining nu |
2140 | | nu.k += 1 |
2141 | | |
2142 | | # "nu[k] := nu[k-1] perp v[k-1]" |
2143 | | nu._clear() |
2144 | | alpha[0] = nu._split_vertex(v[nu.k-1]) |
2145 | | alpha[1] = -1 |
2146 | | i = nu._refine(alpha, n, M, _dig, uif) |
2147 | | if verbosity >= 5: |
2148 | | ST_vis_invariant = int(i) |
2149 | | ST_vis_splitvert = int(v[nu.k-1]) |
2150 | | |
2151 | | mpz_set_si(Lambda_mpz[nu.k], i) |
2152 | | |
2153 | | # only if this is the first time moving down the search tree: |
2154 | | if h == -1: state = 5; continue |
2155 | | |
2156 | | # update hzf |
2157 | | if hzf == nu.k-1 and mpz_cmp(Lambda_mpz[nu.k], zf_mpz[nu.k]) == 0: hzf = nu.k |
2158 | | if not lab: state = 3; continue |
2159 | | |
2160 | | # "qzb := cmp(Lambda[k], zb[k])" |
2161 | | if qzb == 0: |
2162 | | if mpz_cmp_si(zb_mpz[nu.k], -1) == 0: # if "zb[k] == oo" |
2163 | | qzb = -1 |
2164 | | else: |
2165 | | qzb = mpz_cmp( Lambda_mpz[nu.k], zb_mpz[nu.k] ) |
2166 | | # update hzb |
2167 | | if hzb == nu.k-1 and qzb == 0: hzb = nu.k |
2168 | | |
2169 | | # if Lambda[k] > zb[k], then zb[k] := Lambda[k] |
2170 | | # (zb keeps track of the indicator invariants corresponding to |
2171 | | # rho, the closest canonical leaf so far seen- if Lambda is |
2172 | | # bigger, then rho must be about to change |
2173 | | if qzb > 0: mpz_set(zb_mpz[nu.k], Lambda_mpz[nu.k]) |
2174 | | state = 3 |
2175 | | |
2176 | | elif state == 3: # attempt to rule out automorphisms while moving down |
2177 | | # the tree |
2178 | | if hzf <= nu.k or (lab and qzb >= 0): # changed hzb to hzf, == to <= |
2179 | | state = 4 |
2180 | | else: state = 6 |
2181 | | # if k > hzf, then we know that nu currently does not look like |
2182 | | # zeta, the first terminal node encountered. Then if we are not |
2183 | | # looking for a canonical label, there is no reason to continue. |
2184 | | # However, if we are looking for one, and qzb < 0, i.e. |
2185 | | # Lambda[k] < zb[k], then the indicator is not maximal, and we |
2186 | | # can't reach a canonical leaf. |
2187 | | |
2188 | | elif state == 4: # at this point we have -not- ruled out the presence |
2189 | | # of automorphisms |
2190 | | if nu._is_discrete(): state = 7; continue |
2191 | | |
2192 | | # store the first smallest nontrivial cell in W[k], and set v[k] |
2193 | | # equal to its minimum element |
2194 | | v[nu.k] = nu._first_smallest_nontrivial(W[nu.k], n) |
2195 | | |
2196 | | if _dig or not nu._sat_225(n): hh = nu.k + 1 |
2197 | | e[nu.k] = 0 # see state 12, and 17 |
2198 | | state = 2 # continue down the tree |
2199 | | |
2200 | | elif state == 5: # alternative to 3: since we have not yet gotten |
2201 | | # zeta, there are no automorphisms to rule out. |
2202 | | # instead we record Lambda to zf and zb |
2203 | | # (see state 3) |
2204 | | if _base: |
2205 | | base_set.append(v[nu.k-1]) |
2206 | | mpz_set(zf_mpz[nu.k], Lambda_mpz[nu.k]) |
2207 | | mpz_set(zb_mpz[nu.k], Lambda_mpz[nu.k]) |
2208 | | state = 4 |
2209 | | |
2210 | | elif state == 6: # at this stage, there is no reason to continue |
2211 | | # downward, and an automorphism has not been |
2212 | | # discovered |
2213 | | j = nu.k |
2214 | | |
2215 | | # return to the longest ancestor nu[i] of nu that could have a |
2216 | | # descendant equivalent to zeta or could improve on rho. |
2217 | | # All terminal nodes descending from nu[hh] are known to be |
2218 | | # equivalent, so i < hh. Also, if i > hzb, none of the |
2219 | | # descendants of nu[i] can improve rho, since the indicator is |
2220 | | # off (Lambda(nu) < Lambda(rho)). If i >= ht, then no descendant |
2221 | | # of nu[i] is equivalent to zeta (see [1, p67]). |
2222 | | if ht-1 > hzb: |
2223 | | if ht-1 < hh-1: |
2224 | | nu.k = ht-1 |
2225 | | else: |
2226 | | nu.k = hh-1 |
2227 | | else: |
2228 | | if hzb < hh-1: |
2229 | | nu.k = hzb |
2230 | | else: |
2231 | | nu.k = hh-1 |
2232 | | |
2233 | | # TODO: investigate the following line |
2234 | | if nu.k == -1: nu.k = 0 # not in BDM, broke at G = Graph({0:[], 1:[]}), Pi = [[0,1]], lab=False |
2235 | | |
2236 | | if lab: |
2237 | | if hb > nu.k: # update hb since we are backtracking (not in [1]) |
2238 | | hb = nu.k # recall hb is the longest common ancestor of rho and nu |
2239 | | |
2240 | | if j == hh: state = 13; continue |
2241 | | # recall hh: the height of the oldest ancestor of zeta for which |
2242 | | # Lemma 2.25 is satsified, which implies that all terminal nodes |
2243 | | # descended from there are equivalent (or simply k if 2.25 does |
2244 | | # not apply). If we are looking at such a node, then the partition |
2245 | | # at nu[hh] can be used for later pruning, so we store its fixed |
2246 | | # set and a set of representatives of its cells |
2247 | | if l < L-1: l += 1 |
2248 | | for i from 0 <= i < n: |
2249 | | Omega[l][i] = 0 # changed Lambda to Omega |
2250 | | Phi[l][i] = 0 |
2251 | | if nu._is_min_cell_rep(i, hh): |
2252 | | Omega[l][i] = 1 |
2253 | | if nu._is_fixed(i, hh): |
2254 | | Phi[l][i] = 1 |
2255 | | |
2256 | | state = 12 |
2257 | | |
2258 | | elif state == 7: # we have just arrived at a terminal node of the |
2259 | | # search tree T(G, Pi) |
2260 | | # if this is the first terminal node, go directly to 18, to |
2261 | | # process zeta |
2262 | | if h == -1: state = 18; continue |
2263 | | |
2264 | | # hzf is the extremal height of ancestors of both nu and zeta, |
2265 | | # so if k < hzf, nu is not equivalent to zeta, i.e. there is no |
2266 | | # automorphism to discover. |
2267 | | # TODO: investigate why, in practice, the same does not seem to be |
2268 | | # true for hzf < k... BDM had !=, not <, and this broke at |
2269 | | # G = Graph({0:[],1:[],2:[]}), Pi = [[0,1,2]] |
2270 | | if nu.k < hzf: state = 8; continue |
2271 | | |
2272 | | # get the permutation corresponding to this terminal node |
2273 | | nu._get_permutation_from(zeta, gamma) |
2274 | | |
2275 | | if verbosity > 3: |
2276 | | print 'checking for automorphism:' |
2277 | | print [gamma[iii] for iii in range(n)] |
2278 | | |
2279 | | # if G^gamma == G, goto 10 |
2280 | | if _is_automorphism(M, n, gamma): |
2281 | | state = 10 |
2282 | | else: |
2283 | | state = 8 |
2284 | | |
2285 | | elif state == 8: # we have just ruled out the presence of automorphism |
2286 | | # and have not yet considered whether nu improves on |
2287 | | # rho |
2288 | | # if we are not searching for a canonical label, there is nothing |
2289 | | # to do here |
2290 | | if (not lab) or (qzb < 0): |
2291 | | state = 6; continue |
2292 | | |
2293 | | # if Lambda[k] > zb[k] or nu is shorter than rho, then we have |
2294 | | # found an improvement for rho |
2295 | | if (qzb > 0) or (nu.k < k_rho): |
2296 | | state = 9; continue |
2297 | | |
2298 | | # if G(nu) > G(rho) (returns 1), goto 9 |
2299 | | # if G(nu) < G(rho) (returns -1), goto 6 |
2300 | | # if G(nu) == G(rho) (returns 0), get the automorphism and goto 10 |
2301 | | m = nu._compare_with(M, n, rho) |
2302 | | |
2303 | | if m > 0: |
2304 | | state = 9; continue |
2305 | | if m < 0: |
2306 | | state = 6; continue |
2307 | | |
2308 | | rho._get_permutation_from(nu, gamma) |
2309 | | if verbosity > 3: |
2310 | | print 'automorphism discovered:' |
2311 | | print [gamma[iii] for iii in range(n)] |
2312 | | state = 10 |
2313 | | |
2314 | | elif state == 9: # entering this state, nu is a best-so-far guess at |
2315 | | # the canonical label |
2316 | | rho = PartitionStack(nu) |
2317 | | k_rho = nu.k |
2318 | | |
2319 | | qzb = 0 |
2320 | | hb = nu.k |
2321 | | hzb = nu.k |
2322 | | |
2323 | | # set zb[k+1] = Infinity |
2324 | | mpz_set_si(zb_mpz[nu.k+1], -1) |
2325 | | state = 6 |
2326 | | |
2327 | | elif state == 10: # we have an automorphism to process |
2328 | | # increment l |
2329 | | if l < L - 1: |
2330 | | l += 1 |
2331 | | |
2332 | | for i from 0 <= i < n: |
2333 | | if gamma[i] == i: |
2334 | | Phi[l][i] = 1 |
2335 | | Omega[l][i] = 1 |
2336 | | else: |
2337 | | Phi[l][i] = 0 |
2338 | | m = i |
2339 | | j = gamma[i] |
2340 | | while j != i: |
2341 | | if j < m: m = j |
2342 | | j = gamma[j] |
2343 | | if m == i: |
2344 | | Omega[l][i] = 1 |
2345 | | else: |
2346 | | Omega[l][i] = 0 |
2347 | | m = Theta._incorporate_permutation(gamma, n) |
2348 | | # if each orbit of gamma is part of an orbit in Theta, then the |
2349 | | # automorphism is already in the span of those we have seen |
2350 | | if not m: |
2351 | | state = 11 |
2352 | | continue |
2353 | | |
2354 | | # record the automorphism |
2355 | | output.append([ Integer(gamma[i]) for i from 0 <= i < n ]) |
2356 | | |
2357 | | # The variable tvh represents the minimum element of W[k], |
2358 | | # the last time we were at state 13 and backtracking up |
2359 | | # zeta. If this is not still a minimal cell representative of Theta, |
2360 | | # then we need to immediately backtrack to the place where it was |
2361 | | # defined on a part of zeta, since the rest of the tree is now |
2362 | | # equivalent. Otherwise, proceed to 11 and 12 before moving back to |
2363 | | # the hub state. |
2364 | | if Theta.elements[tvh] == -1: |
2365 | | state = 11 |
2366 | | continue |
2367 | | nu.k = h |
2368 | | state = 13 |
2369 | | |
2370 | | elif state == 11: # we have just discovered an automorphism, |
2371 | | # but tvh is still a minimal representative for its |
2372 | | # orbit in Theta. Therefore we cannot backtrack all |
2373 | | # the way to where zeta meets nu. Instead we just use |
2374 | | # indicator values to determine where to backtrack. |
2375 | | if lab: |
2376 | | nu.k = hb |
2377 | | else: |
2378 | | nu.k = h |
2379 | | state = 12 |
2380 | | |
2381 | | elif state == 12: |
2382 | | # e keeps track of the whether W[k] has been thinned out by Omega |
2383 | | # and Phi. It is set to 1 when you have just finished coming up the |
2384 | | # search tree, and have intersected W[k] with Omega[i], for the |
2385 | | # appropriate i < l, but since there may be an automorphism mapping |
2386 | | # one element of W[k] to another still, we thin out W[k] again. |
2387 | | # (see state 17) Coming from 11, this is an explicit automorphism. |
2388 | | # Coming from 6, this is an implicit automorphism. |
2389 | | if e[nu.k] == 1: |
2390 | | for j from 0 <= j < n: |
2391 | | if W[nu.k][j] and not Omega[l][j]: |
2392 | | W[nu.k][j] = 0 |
2393 | | state = 13 |
2394 | | |
2395 | | elif state == 13: # hub state |
2396 | | if nu.k == -1: |
2397 | | # the algorithm has finished |
2398 | | state = -1; continue |
2399 | | if nu.k > h: |
2400 | | # if we are not at a node of zeta |
2401 | | state = 17; continue |
2402 | | if nu.k == h: |
2403 | | # if we are at a node of zeta, then we have not yet backtracked |
2404 | | # UP zeta, so skip the rest of 13 |
2405 | | state = 14; continue |
2406 | | |
2407 | | # thus, it must be that k < h, and this means we are done |
2408 | | # searching underneath zeta[k+1], so now, k is the new longest |
2409 | | # ancestor of nu and zeta: |
2410 | | h = nu.k |
2411 | | |
2412 | | # set tvh to the minimum cell representative of W[k] |
2413 | | # (see states 10 and 14) |
2414 | | for i from 0 <= i < n: |
2415 | | if W[nu.k][i]: |
2416 | | tvh = i |
2417 | | break |
2418 | | state = 14 |
2419 | | |
2420 | | elif state == 14: # iterate v[k] through W[k] until a minimum cell rep |
2421 | | # of Theta is found: |
2422 | | # this state gets hit only when we are looking for a |
2423 | | # new split off of zeta |
2424 | | # The variable tvh was set to be the minimum element of W[k] |
2425 | | # the last time we were at state 13 and backtracking up |
2426 | | # zeta. If this is in the same cell of Theta as v[k], increment |
2427 | | # index (see Theorem 2.33 in [1]) |
2428 | | if Theta._find(v[nu.k]) == Theta._find(tvh): |
2429 | | index += 1 |
2430 | | |
2431 | | # find the next v[k] in W[k] |
2432 | | i = v[nu.k] + 1 |
2433 | | while i < n and not W[nu.k][i]: |
2434 | | i += 1 |
2435 | | if i < n: |
2436 | | v[nu.k] = i |
2437 | | else: |
2438 | | # there is no new vertex to consider at this level |
2439 | | v[nu.k] = -1 |
2440 | | state = 16 |
2441 | | continue |
2442 | | |
2443 | | # if the new v[k] is not a minimum cell representative of Theta, |
2444 | | # then we already considered that rep., and that subtree was |
2445 | | # isomorphic to the one corresponding to v[k] |
2446 | | if Theta.elements[v[nu.k]] != -1: state = 14 |
2447 | | else: |
2448 | | # otherwise, we do have a vertex to consider |
2449 | | state = 15 |
2450 | | |
2451 | | elif state == 15: # we have a new vertex, v[k], that we must split on |
2452 | | # hh is smallest such that nu[hh] satisfies Lemma 2.25. If it is |
2453 | | # larger than k+1, it must be modified, since we are changing that |
2454 | | # part |
2455 | | if nu.k + 1 < hh: |
2456 | | hh = nu.k + 1 |
2457 | | # hzf is maximal such that indicators line up for nu and zeta |
2458 | | if nu.k < hzf: |
2459 | | hzf = nu.k |
2460 | | if not lab or hzb < nu.k: |
2461 | | # in either case there is no need to update hzb, which is the |
2462 | | # length for which nu and rho have the same indicators |
2463 | | state = 2; continue |
2464 | | hzb = nu.k |
2465 | | qzb = 0 |
2466 | | state = 2 |
2467 | | |
2468 | | elif state == 16: # backtrack one level in the search tree, recording |
2469 | | # information relevant to Theorem 2.33 |
2470 | | j = 0 |
2471 | | for i from 0 <= i < n: |
2472 | | if W[nu.k][i]: j += 1 |
2473 | | if j == index and ht == nu.k+1: ht = nu.k |
2474 | | size *= index |
2475 | | index = 0 |
2476 | | nu.k -= 1 |
2477 | | |
2478 | | if lab: |
2479 | | if hb > nu.k: # update hb since we are backtracking (not in [1]): |
2480 | | hb = nu.k # recall hb is the longest common ancestor of rho and nu |
2481 | | |
2482 | | state = 13 |
2483 | | |
2484 | | elif state == 17: # you have just finished coming up the search tree, |
2485 | | # and must now consider going back down. |
2486 | | if e[nu.k] == 0: |
2487 | | # intersect W[k] with each Omega[i] such that {v_0,...,v_(k-1)} |
2488 | | # is contained in Phi[i] |
2489 | | for i from 0 <= i <= l: |
2490 | | # check if {v_0,...,v_(k-1)} is contained in Phi[i] |
2491 | | # i.e. fixed pointwise by the automorphisms so far seen |
2492 | | j = 0 |
2493 | | while j < nu.k and Phi[i][v[j]]: |
2494 | | j += 1 |
2495 | | # if so, only check the minimal orbit representatives |
2496 | | if j == nu.k: |
2497 | | for j from 0 <= j < n: |
2498 | | if W[nu.k][j] and not Omega[i][j]: |
2499 | | W[nu.k][j] = 0 |
2500 | | e[nu.k] = 1 # see state 12 |
2501 | | |
2502 | | # see if there is a relevant vertex to split on: |
2503 | | i = v[nu.k] + 1 |
2504 | | while i < n and not W[nu.k][i]: |
2505 | | i += 1 |
2506 | | if i < n: |
2507 | | v[nu.k] = i |
2508 | | state = 15 |
2509 | | continue |
2510 | | else: |
2511 | | v[nu.k] = -1 |
2512 | | |
2513 | | # otherwise backtrack one level |
2514 | | nu.k -= 1 |
2515 | | state = 13 |
2516 | | |
2517 | | elif state == 18: # The first time we encounter a terminal node, we |
2518 | | # come straight here to set up zeta. This is a one- |
2519 | | # time state. |
2520 | | # initialize counters for zeta: |
2521 | | h = nu.k # zeta[h] == nu[h] |
2522 | | ht = nu.k # nodes descended from zeta[ht] are all equivalent |
2523 | | hzf = nu.k # max such that indicators for zeta and nu agree |
2524 | | |
2525 | | zeta = PartitionStack(nu) |
2526 | | |
2527 | | nu.k -= 1 |
2528 | | if not lab: state = 13; continue |
2529 | | |
2530 | | rho = PartitionStack(nu) |
2531 | | |
2532 | | # initialize counters for rho: |
2533 | | k_rho = nu.k + 1 # number of partitions in rho |
2534 | | hzb = nu.k # max such that indicators for rho and nu agree - BDM had k+1 |
2535 | | hb = nu.k # rho[hb] == nu[hb] - BDM had k+1 |
2536 | | |
2537 | | qzb = 0 # Lambda[k] == zb[k], so... |
2538 | | state = 13 |
2539 | | |
2540 | | # free the GMP ints |
2541 | | for i from 0 <= i < n+2: |
2542 | | mpz_clear(Lambda_mpz[i]) |
2543 | | mpz_clear(zf_mpz[i]) |
2544 | | mpz_clear(zb_mpz[i]) |
2545 | | |
2546 | | # free int arrays |
2547 | | sage_free(gamma) |
2548 | | sage_free(W[0]) |
2549 | | sage_free(Phi[0]) |
2550 | | sage_free(Omega[0]) |
2551 | | sage_free(alpha) |
2552 | | sage_free(v) |
2553 | | sage_free(e) |
2554 | | |
2555 | | # free GMP int pointers |
2556 | | sage_free(Lambda_mpz) |
2557 | | sage_free(zf_mpz) |
2558 | | sage_free(zb_mpz) |
2559 | | |
2560 | | # free int pointers |
2561 | | sage_free(W) |
2562 | | sage_free(Phi) |
2563 | | sage_free(Omega) |
2564 | | |
2565 | | # prepare output |
2566 | | if not (lab or dict_rep or certify or base or order): |
2567 | | return output |
2568 | | output_tuple = [output] |
2569 | | if dict_rep: |
2570 | | if isinstance(G, GenericGraph): |
2571 | | ddd = {} |
2572 | | for vvv in ffrom.iterkeys(): # v is a C variable |
2573 | | if ffrom[vvv] != 0: |
2574 | | ddd[vvv] = ffrom[vvv] |
2575 | | else: |
2576 | | ddd[vvv] = n |
2577 | | else: # G is a CGraph |
2578 | | ddd = {} |
2579 | | for i from 0 <= i < n: |
2580 | | ddd[i] = i |
2581 | | output_tuple.append(ddd) |
2582 | | if lab: |
2583 | | H = _term_pnest_graph(G, rho) |
2584 | | output_tuple.append(H) |
2585 | | if certify: |
2586 | | dd = {} |
2587 | | for i from 0 <= i < n: |
2588 | | dd[to[rho.entries[i]]] = i |
2589 | | output_tuple.append(dd) |
2590 | | if base: |
2591 | | output_tuple.append(base_set) |
2592 | | if order: |
2593 | | output_tuple.append(size) |
2594 | | return tuple(output_tuple) |
2595 | | |
2596 | | # Benchmarking functions |
2597 | | |
2598 | | def all_labeled_graphs(n): |
2599 | | """ |
2600 | | Returns all labeled graphs on n vertices 0,1,...,n-1. Used in |
2601 | | classifying isomorphism types (naive approach), and more |
2602 | | importantly in benchmarking the search algorithm. |
2603 | | |
2604 | | EXAMPLE:: |
2605 | | |
2606 | | sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree |
2607 | | sage: from sage.graphs.graph_isom import all_labeled_graphs |
2608 | | sage: Glist = {} |
2609 | | sage: Giso = {} |
2610 | | sage: for n in range(1,5): |
2611 | | ... Glist[n] = all_labeled_graphs(n) |
2612 | | ... Giso[n] = [] |
2613 | | ... for g in Glist[n]: |
2614 | | ... a, b = search_tree(g, [range(n)]) |
2615 | | ... inn = False |
2616 | | ... for gi in Giso[n]: |
2617 | | ... if b == gi: |
2618 | | ... inn = True |
2619 | | ... if not inn: |
2620 | | ... Giso[n].append(b) |
2621 | | sage: for n in Giso: |
2622 | | ... print n, len(Giso[n]) |
2623 | | 1 1 |
2624 | | 2 2 |
2625 | | 3 4 |
2626 | | 4 11 |
2627 | | sage: n = 5 |
2628 | | sage: Glist[n] = all_labeled_graphs(n) |
2629 | | sage: Giso[n] = [] |
2630 | | sage: for g in Glist[5]: |
2631 | | ... a, b = search_tree(g, [range(n)]) |
2632 | | ... inn = False |
2633 | | ... for gi in Giso[n]: |
2634 | | ... if b == gi: |
2635 | | ... inn = True |
2636 | | ... if not inn: |
2637 | | ... Giso[n].append(b) |
2638 | | sage: print n, len(Giso[n]) # long time |
2639 | | 5 34 |
2640 | | sage: graphs_list.show_graphs(Giso[4]) |
2641 | | """ |
2642 | | TE = [] |
2643 | | for i in range(n): |
2644 | | for j in range(i): |
2645 | | TE.append((i, j)) |
2646 | | m = len(TE) |
2647 | | Glist= [] |
2648 | | for i in range(2**m): |
2649 | | G = Graph(n) |
2650 | | b = Integer(i).binary() |
2651 | | b = '0'*(m-len(b)) + b |
2652 | | for i in range(m): |
2653 | | if int(b[i]): |
2654 | | G.add_edge(TE[i]) |
2655 | | Glist.append(G) |
2656 | | return Glist |
2657 | | |
2658 | | def kpow(listy, k): |
2659 | | """ |
2660 | | Returns the subset of the power set of listy consisting of subsets |
2661 | | of size k. Used in all_ordered_partitions. |
2662 | | |
2663 | | EXAMPLE:: |
2664 | | |
2665 | | sage: from sage.graphs.graph_isom import kpow |
2666 | | sage: kpow(['a', 1, {}], 2) |
2667 | | [[1, 'a'], [{}, 'a'], ['a', 1], [{}, 1], ['a', {}], [1, {}]] |
2668 | | """ |
2669 | | list = [] |
2670 | | if k > 1: |
2671 | | for LL in kpow(listy, k-1): |
2672 | | for a in listy: |
2673 | | if not a in LL: |
2674 | | list.append([a] + LL) |
2675 | | if k == 1: |
2676 | | for i in listy: |
2677 | | list.append([i]) |
2678 | | return list |
2679 | | |
2680 | | def all_ordered_partitions(listy): |
2681 | | """ |
2682 | | Returns all ordered partitions of the set 0,1,...,n-1. Used in |
2683 | | benchmarking the search algorithm. |
2684 | | |
2685 | | EXAMPLE:: |
2686 | | |
2687 | | sage: from sage.graphs.graph_isom import all_ordered_partitions |
2688 | | sage: all_ordered_partitions(['a', 1, {}]) |
2689 | | [[['a'], [1], [{}]], |
2690 | | [['a'], [{}], [1]], |
2691 | | [['a'], [{}, 1]], |
2692 | | [['a'], [1, {}]], |
2693 | | [[1], ['a'], [{}]], |
2694 | | [[1], [{}], ['a']], |
2695 | | [[1], [{}, 'a']], |
2696 | | [[1], ['a', {}]], |
2697 | | [[{}], ['a'], [1]], |
2698 | | [[{}], [1], ['a']], |
2699 | | [[{}], [1, 'a']], |
2700 | | [[{}], ['a', 1]], |
2701 | | [[1, 'a'], [{}]], |
2702 | | [[{}, 'a'], [1]], |
2703 | | [['a', 1], [{}]], |
2704 | | [[{}, 1], ['a']], |
2705 | | [['a', {}], [1]], |
2706 | | [[1, {}], ['a']], |
2707 | | [[{}, 1, 'a']], |
2708 | | [[1, {}, 'a']], |
2709 | | [[{}, 'a', 1]], |
2710 | | [['a', {}, 1]], |
2711 | | [[1, 'a', {}]], |
2712 | | [['a', 1, {}]]] |
2713 | | """ |
2714 | | LL = [] |
2715 | | for i in range(1,len(listy)+1): |
2716 | | for cell in kpow(listy, i): |
2717 | | list_remainder = [x for x in listy if x not in cell] |
2718 | | remainder_partitions = all_ordered_partitions(list_remainder) |
2719 | | for remainder in remainder_partitions: |
2720 | | LL.append( [cell] + remainder ) |
2721 | | if len(listy) == 0: |
2722 | | return [[]] |
2723 | | else: |
2724 | | return LL |
2725 | | |
2726 | | def all_labeled_digraphs_with_loops(n): |
2727 | | """ |
2728 | | Returns all labeled digraphs on n vertices 0,1,...,n-1. Used in |
2729 | | classifying isomorphism types (naive approach), and more |
2730 | | importantly in benchmarking the search algorithm. |
2731 | | |
2732 | | EXAMPLE:: |
2733 | | |
2734 | | sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree |
2735 | | sage: from sage.graphs.graph_isom import all_labeled_digraphs_with_loops |
2736 | | sage: Glist = {} |
2737 | | sage: Giso = {} |
2738 | | sage: for n in range(1,4): |
2739 | | ... Glist[n] = all_labeled_digraphs_with_loops(n) |
2740 | | ... Giso[n] = [] |
2741 | | ... for g in Glist[n]: |
2742 | | ... a, b = search_tree(g, [range(n)], dig=True) |
2743 | | ... inn = False |
2744 | | ... for gi in Giso[n]: |
2745 | | ... if b == gi: |
2746 | | ... inn = True |
2747 | | ... if not inn: |
2748 | | ... Giso[n].append(b) |
2749 | | sage: for n in Giso: |
2750 | | ... print n, len(Giso[n]) |
2751 | | 1 2 |
2752 | | 2 10 |
2753 | | 3 104 |
2754 | | """ |
2755 | | TE = [] |
2756 | | for i in range(n): |
2757 | | for j in range(n): |
2758 | | TE.append((i, j)) |
2759 | | m = len(TE) |
2760 | | Glist= [] |
2761 | | for i in range(2**m): |
2762 | | G = DiGraph(n, loops=True) |
2763 | | b = Integer(i).binary() |
2764 | | b = '0'*(m-len(b)) + b |
2765 | | for j in range(m): |
2766 | | if int(b[j]): |
2767 | | G.add_edge(TE[j]) |
2768 | | Glist.append(G) |
2769 | | return Glist |
2770 | | |
2771 | | def all_labeled_digraphs(n): |
2772 | | """ |
2773 | | EXAMPLES:: |
2774 | | |
2775 | | sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree |
2776 | | sage: from sage.graphs.graph_isom import all_labeled_digraphs |
2777 | | sage: Glist = {} |
2778 | | sage: Giso = {} |
2779 | | sage: for n in range(1,4): |
2780 | | ... Glist[n] = all_labeled_digraphs(n) |
2781 | | ... Giso[n] = [] |
2782 | | ... for g in Glist[n]: |
2783 | | ... a, b = search_tree(g, [range(n)], dig=True) |
2784 | | ... inn = False |
2785 | | ... for gi in Giso[n]: |
2786 | | ... if b == gi: |
2787 | | ... inn = True |
2788 | | ... if not inn: |
2789 | | ... Giso[n].append(b) |
2790 | | sage: for n in Giso: |
2791 | | ... print n, len(Giso[n]) |
2792 | | 1 1 |
2793 | | 2 3 |
2794 | | 3 16 |
2795 | | """ |
2796 | | |
2797 | | ## sage: n = 4 # long time (4 minutes) |
2798 | | ## sage: Glist[n] = all_labeled_digraphs(n) # long time |
2799 | | ## sage: Giso[n] = [] # long time |
2800 | | ## sage: for g in Glist[n]: # long time |
2801 | | ## ... a, b = search_tree(g, [range(n)], dig=True) |
2802 | | ## ... inn = False |
2803 | | ## ... for gi in Giso[n]: |
2804 | | ## ... if enum(b) == enum(gi): |
2805 | | ## ... inn = True |
2806 | | ## ... if not inn: |
2807 | | ## ... Giso[n].append(b) |
2808 | | ## sage: print n, len(Giso[n]) # long time |
2809 | | ## 4 218 |
2810 | | TE = [] |
2811 | | for i in range(n): |
2812 | | for j in range(n): |
2813 | | if i != j: TE.append((i, j)) |
2814 | | m = len(TE) |
2815 | | Glist= [] |
2816 | | for i in range(2**m): |
2817 | | G = DiGraph(n, loops=True) |
2818 | | b = Integer(i).binary() |
2819 | | b = '0'*(m-len(b)) + b |
2820 | | for j in range(m): |
2821 | | if int(b[j]): |
2822 | | G.add_edge(TE[j]) |
2823 | | Glist.append(G) |
2824 | | return Glist |
2825 | | |
2826 | | def perm_group_elt(lperm): |
2827 | | """ |
2828 | | Given a list permutation of the set 0, 1, ..., n-1, returns the |
2829 | | corresponding PermutationGroupElement where we take 0 = n. |
2830 | | |
2831 | | EXAMPLE:: |
2832 | | |
2833 | | sage: from sage.graphs.graph_isom import perm_group_elt |
2834 | | sage: perm_group_elt([0,2,1]) |
2835 | | (1,2) |
2836 | | sage: perm_group_elt([1,2,0]) |
2837 | | (1,2,3) |
2838 | | """ |
2839 | | from sage.groups.perm_gps.permgroup_named import SymmetricGroup |
2840 | | n = len(lperm) |
2841 | | S = SymmetricGroup(n) |
2842 | | Part = orbit_partition(lperm, list_perm=True) |
2843 | | gens = [] |
2844 | | for z in Part: |
2845 | | if len(z) > 1: |
2846 | | if 0 in z: |
2847 | | zed = z.index(0) |
2848 | | generator = z[:zed] + [n] + z[zed+1:] |
2849 | | gens.append(tuple(generator)) |
2850 | | else: |
2851 | | gens.append(tuple(z)) |
2852 | | E = S(gens) |
2853 | | return E |
2854 | | |
2855 | | def orbit_partition(gamma, list_perm=False): |
2856 | | r""" |
2857 | | Assuming that G is a graph on vertices 0,1,...,n-1, and gamma is an |
2858 | | element of SymmetricGroup(n), returns the partition of the vertex |
2859 | | set determined by the orbits of gamma, considered as action on the |
2860 | | set 1,2,...,n where we take 0 = n. In other words, returns the |
2861 | | partition determined by a cyclic representation of gamma. |
2862 | | |
2863 | | INPUT: |
2864 | | |
2865 | | |
2866 | | - ``list_perm`` - if True, assumes |
2867 | | ``gamma`` is a list representing the map |
2868 | | `i \mapsto ``gamma``[i]`. |
2869 | | |
2870 | | |
2871 | | EXAMPLES:: |
2872 | | |
2873 | | sage: from sage.graphs.graph_isom import orbit_partition |
2874 | | sage: G = graphs.PetersenGraph() |
2875 | | sage: S = SymmetricGroup(10) |
2876 | | sage: gamma = S('(10,1,2,3,4)(5,6,7)(8,9)') |
2877 | | sage: orbit_partition(gamma) |
2878 | | [[1, 2, 3, 4, 0], [5, 6, 7], [8, 9]] |
2879 | | sage: gamma = S('(10,5)(1,6)(2,7)(3,8)(4,9)') |
2880 | | sage: orbit_partition(gamma) |
2881 | | [[1, 6], [2, 7], [3, 8], [4, 9], [5, 0]] |
2882 | | """ |
2883 | | if list_perm: |
2884 | | n = len(gamma) |
2885 | | seen = [1] + [0]*(n-1) |
2886 | | i = 0 |
2887 | | p = 0 |
2888 | | partition = [[0]] |
2889 | | while sum(seen) < n: |
2890 | | if gamma[i] != partition[p][0]: |
2891 | | partition[p].append(gamma[i]) |
2892 | | i = gamma[i] |
2893 | | seen[i] = 1 |
2894 | | else: |
2895 | | i = min([j for j in range(n) if seen[j] == 0]) |
2896 | | partition.append([i]) |
2897 | | p += 1 |
2898 | | seen[i] = 1 |
2899 | | return partition |
2900 | | else: |
2901 | | n = len(gamma.list()) |
2902 | | l = [] |
2903 | | for i in range(1,n+1): |
2904 | | orb = gamma.orbit(i) |
2905 | | if orb not in l: l.append(orb) |
2906 | | for i in l: |
2907 | | for j in range(len(i)): |
2908 | | if i[j] == n: |
2909 | | i[j] = 0 |
2910 | | return l |
2911 | | |
2912 | | def verify_partition_refinement(G, initial_partition, terminal_partition): |
2913 | | """ |
2914 | | Verify that the refinement is correct. |
2915 | | |
2916 | | EXAMPLE:: |
2917 | | |
2918 | | sage: from sage.graphs.graph_isom import PartitionStack |
2919 | | sage: from sage.graphs.base.sparse_graph import SparseGraph |
2920 | | sage: G = SparseGraph(10) |
2921 | | sage: for i,j,_ in graphs.PetersenGraph().edge_iterator(): |
2922 | | ... G.add_arc(i,j) |
2923 | | ... G.add_arc(j,i) |
2924 | | sage: P = PartitionStack(10) |
2925 | | sage: P.set_k(1) |
2926 | | sage: P.split_vertex(0) |
2927 | | sage: P.refine(G, [0], 10, 0, 1) |
2928 | | sage: P |
2929 | | (0,2,3,6,7,8,9,1,4,5) |
2930 | | (0|2,3,6,7,8,9|1,4,5) |
2931 | | sage: P.set_k(2) |
2932 | | sage: P.split_vertex(1) |
2933 | | |
2934 | | Note that this line implicitly tests the function |
2935 | | verify_partition_refinement:: |
2936 | | |
2937 | | sage: P.refine(G, [7], 10, 0, 1, test=True) |
2938 | | sage: P |
2939 | | (0,3,7,8,9,2,6,1,4,5) |
2940 | | (0|3,7,8,9,2,6|1,4,5) |
2941 | | (0|3,7,8,9|2,6|1|4,5) |
2942 | | """ |
2943 | | if not G.is_equitable(terminal_partition): |
2944 | | raise RuntimeError("Resulting partition is not equitable!!!!!!!!!\n"+\ |
2945 | | str(initial_partition) + "\n" + \ |
2946 | | str(terminal_partition) + "\n" + \ |
2947 | | str(G.am())) |
2948 | | |