sage-main
view sage/graphs/graph_isom.pyx @ 11802:5c72dbb92d82
3.4
| author | mabshoff@sage.math.washington.edu |
|---|---|
| date | Wed Mar 11 21:57:15 2009 -0700 (17 months ago) |
| parents | 093fe93e7563 |
| children | a22d9145c761 |
line source
1 """
2 N.I.C.E. - Nice (as in open source) Isomorphism Check Engine
4 Automorphism group computation and isomorphism checking for
5 graphs.
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.
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.
35 AUTHORS:
37 - Robert L. Miller (2007-03-20): initial version
39 - Tom Boothby (2007-03-20): help with indicator function
41 - Robert L. Miller (2007-04-07-30): optimizations
43 - Robert L. Miller (2007-07-07-14): PartitionStack and OrbitPartition
45 - Tom Boothby (2007-07-14) datastructure advice
47 - Robert L. Miller (2007-07-16-20): bug fixes
49 REFERENCE:
51 - [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus
52 Numerantium, Vol. 30 (1981), pp. 45-87.
54 .. note::
56 #. Often we assume that G is a graph on vertices
57 0,1,...,n-1.
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 """
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 #*****************************************************************************
71 from sage.graphs.graph import GenericGraph, Graph, DiGraph
72 from sage.misc.misc import cputime
73 from sage.rings.integer cimport Integer
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.
81 EXAMPLES::
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
96 ::
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 """
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
134 def __dealloc__(self):
135 sage_free(self.elements)
136 sage_free(self.sizes)
138 def find(self, x):
139 """
140 Returns an element of the cell which depends only on the cell.
142 EXAMPLE::
144 sage: from sage.graphs.graph_isom import OrbitPartition
145 sage: K = OrbitPartition(20)
147 0 and 1 begin in different cells::
149 sage: K.find(0)
150 0
151 sage: K.find(1)
152 1
154 Now we put them in the same cell::
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)
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]
170 def union_find(self, a, b):
171 """
172 Merges the cells containing a and b.
174 EXAMPLE::
176 sage: from sage.graphs.graph_isom import OrbitPartition
177 sage: K = OrbitPartition(20)
179 0 and 1 begin in different cells::
181 sage: K.find(0)
182 0
183 sage: K.find(1)
184 1
186 Now we put them in the same cell::
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)
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
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]
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.
216 EXAMPLE::
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
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
240 def vee_with(self, other, n):
241 """
242 Merges the minimal number of cells such that other is finer than
243 self.
245 EXAMPLE::
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)
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)))
267 cdef int _is_min_cell_rep(self, int i):
268 if self.elements[i] == -1:
269 return 1
270 return 0
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
277 def incorporate_permutation(self, gamma):
278 """
279 Unions the cells of self which contain common elements of some
280 orbit of gamma.
282 INPUT:
285 - ``gamma`` - a permutation, in list notation
288 EXAMPLE::
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)
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
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
331 cdef class PartitionStack:
332 """
333 TODO: documentation
335 EXAMPLES::
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
376 ::
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.")
456 def __dealloc__(self):
457 sage_free(self.entries)
458 sage_free(self.levels)
460 def __repr__(self):
461 """
462 EXAMPLE::
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
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.
509 EXAMPLE::
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)
544 ::
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
563 def set_k(self, k):
564 """
565 Sets self.k, the index of the finest partition.
567 ::
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)
602 ::
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
612 def is_discrete(self):
613 """
614 Returns whether the partition consists of only singletons.
616 EXAMPLE::
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)
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
671 def num_cells(self):
672 """
673 Return the number of cells in the finest partition.
675 EXAMPLE::
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()
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
708 def sat_225(self, n):
709 """
710 Whether the finest partition satisfies the hypotheses of Lemma 2.25
711 in [1].
713 EXAMPLE::
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
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
772 cdef int _is_min_cell_rep(self, int i, int k):
773 return i == 0 or self.levels[i-1] <= k
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
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).
786 EXAMPLE::
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)
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
815 def percolate(self, start, end):
816 """
817 Perform one round of bubble sort, moving the smallest element to
818 the front.
820 EXAMPLE::
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)
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
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.
847 EXAMPLE::
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
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)])
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
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]
897 for j from i >= j >= 0:
898 counts[degrees[j]] -= 1
899 output[counts[degrees[j]]] = self.entries[start+j]
901 max_location = counts[max_location]+start
903 for j from 0 <= j <= i:
904 self.entries[start+j] = output[j]
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
913 return max_location
915 def clear(self):
916 """
917 Merges all cells in the partition stack.
919 EXAMPLE::
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()
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
952 def refine(self, CGraph G, alpha, n, dig, uif, test=False):
953 """
954 Implementation of Algorithm 2.5 in [1].
956 EXAMPLE::
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
997 ::
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)
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
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
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.
1192 EXAMPLE::
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
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
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
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
1259 cdef void _get_permutation_from(self, PartitionStack zeta, int *gamma):
1260 cdef int i = 0
1262 while True:
1263 gamma[zeta.entries[i]] = self.entries[i]
1264 i += 1
1265 if self.levels[i-1] == -1: break
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
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
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).
1293 EXAMPLE::
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
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.
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.
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
1377 - ``use_indicator_function`` - option to turn off
1378 indicator function (False - slower)
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)
1384 - ``base`` - whether to return the first sequence of
1385 split vertices (used in computing the order of the group)
1387 - ``order`` - whether to return the order of the
1388 automorphism group
1391 STATE DIAGRAM::
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)
1409 .. note::
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).
1416 EXAMPLES: The following example is due to Chris Godsil::
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
1428 ::
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
1436 ::
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 + 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 3]
1474 [ 3 -1 0 -1 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 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 sage: DodecAut2 = PermutationGroup([perm_group_elt(cc) for cc in c])
1482 sage: DodecAut2.character_table()
1483 [ 1 1 1 1 1 1 1 1 1 1]
1484 [ 1 -1 1 1 -1 1 -1 1 -1 -1]
1485 [ 3 -1 0 -1 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 3]
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 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 -zeta5^3 - zeta5^2 - 1 -3]
1488 [ 3 1 0 -1 -zeta5^3 - zeta5^2 - 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 zeta5^3 + zeta5^2 -3]
1489 [ 4 0 1 0 -1 -1 1 -1 -1 4]
1490 [ 4 0 1 0 1 -1 -1 -1 1 -4]
1491 [ 5 1 -1 1 0 0 -1 0 0 5]
1492 [ 5 -1 -1 1 0 0 1 0 0 -5]
1494 ::
1496 sage: G = graphs.PetersenGraph()
1497 sage: Pi=[range(10)]
1498 sage: a,b = search_tree(G, Pi)
1499 sage: print a, b.graph6_string()
1500 [[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?
1501 sage: c = search_tree(G, Pi, lab=False)
1502 sage: PAut = PermutationGroup([perm_group_elt(aa) for aa in a])
1503 sage: PAut.character_table()
1504 [ 1 1 1 1 1 1 1]
1505 [ 1 -1 1 -1 1 -1 1]
1506 [ 4 -2 0 1 1 0 -1]
1507 [ 4 2 0 -1 1 0 -1]
1508 [ 5 1 1 1 -1 -1 0]
1509 [ 5 -1 1 -1 -1 1 0]
1510 [ 6 0 -2 0 0 0 1]
1511 sage: PAut = PermutationGroup([perm_group_elt(cc) for cc in c])
1512 sage: PAut.character_table()
1513 [ 1 1 1 1 1 1 1]
1514 [ 1 -1 1 -1 1 -1 1]
1515 [ 4 -2 0 1 1 0 -1]
1516 [ 4 2 0 -1 1 0 -1]
1517 [ 5 1 1 1 -1 -1 0]
1518 [ 5 -1 1 -1 -1 1 0]
1519 [ 6 0 -2 0 0 0 1]
1521 ::
1523 sage: G = graphs.CubeGraph(3)
1524 sage: Pi = []
1525 sage: for i in range(8):
1526 ... b = Integer(i).binary()
1527 ... Pi.append(b.zfill(3))
1528 ...
1529 sage: Pi = [Pi]
1530 sage: a,b = search_tree(G, Pi)
1531 sage: print a, b.graph6_string()
1532 [[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_
1533 sage: c = search_tree(G, Pi, lab=False)
1535 ::
1537 sage: PermutationGroup([perm_group_elt(aa) for aa in a]).order()
1538 48
1539 sage: PermutationGroup([perm_group_elt(cc) for cc in c]).order()
1540 48
1541 sage: DodecAut.order()
1542 120
1543 sage: PAut.order()
1544 120
1546 ::
1548 sage: D = graphs.DodecahedralGraph()
1549 sage: a,b,c = search_tree(D, [range(20)], certify=True)
1550 sage: from sage.plot.plot import GraphicsArray
1551 sage: from sage.graphs.graph_fast import spring_layout_fast
1552 sage: position_D = spring_layout_fast(D)
1553 sage: position_b = {}
1554 sage: for vert in position_D:
1555 ... position_b[c[vert]] = position_D[vert]
1556 sage: graphics_array([D.plot(pos=position_D), b.plot(pos=position_b)]).show()
1557 sage: c
1558 {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}
1560 BENCHMARKS: The following examples are given to check modifications
1561 to the algorithm for optimization.
1563 ::
1565 sage: G = Graph({0:[]})
1566 sage: Pi = [[0]]
1567 sage: a,b = search_tree(G, Pi)
1568 sage: print a, b.graph6_string()
1569 [] @
1570 sage: a,b = search_tree(G, Pi, dig=True)
1571 sage: print a, b.graph6_string()
1572 [] @
1573 sage: search_tree(G, Pi, lab=False)
1574 []
1576 ::
1578 sage: from sage.graphs.graph_isom import all_labeled_graphs, all_ordered_partitions
1580 ::
1582 sage: graph2 = all_labeled_graphs(2)
1583 sage: part2 = all_ordered_partitions(range(2))
1584 sage: for G in graph2:
1585 ... for Pi in part2:
1586 ... a,b = search_tree(G, Pi)
1587 ... c,d = search_tree(G, Pi, dig=True)
1588 ... e = search_tree(G, Pi, lab=False)
1589 ... a = str(a); b = b.graph6_string(); c = str(c); d = d.graph6_string(); e = str(e)
1590 ... print a.ljust(15), b.ljust(5), c.ljust(15), d.ljust(5), e.ljust(15)
1591 [] A? [] A? []
1592 [] A? [] A? []
1593 [[1, 0]] A? [[1, 0]] A? [[1, 0]]
1594 [[1, 0]] A? [[1, 0]] A? [[1, 0]]
1595 [] A_ [] A_ []
1596 [] A_ [] A_ []
1597 [[1, 0]] A_ [[1, 0]] A_ [[1, 0]]
1598 [[1, 0]] A_ [[1, 0]] A_ [[1, 0]]
1600 ::
1602 sage: graph3 = all_labeled_graphs(3)
1603 sage: part3 = all_ordered_partitions(range(3))
1604 sage: for G in graph3:
1605 ... for Pi in part3:
1606 ... a,b = search_tree(G, Pi)
1607 ... c,d = search_tree(G, Pi, dig=True)
1608 ... e = search_tree(G, Pi, lab=False)
1609 ... a = str(a); b = b.graph6_string(); c = str(c); d = d.graph6_string(); e = str(e)
1610 ... print a.ljust(15), b.ljust(5), c.ljust(15), d.ljust(5), e.ljust(15)
1611 [] B? [] B? []
1612 [] B? [] B? []
1613 [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]]
1614 [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]]
1615 [] B? [] B? []
1616 [] B? [] B? []
1617 [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]]
1618 [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]]
1619 [] B? [] B? []
1620 [] B? [] B? []
1621 [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]]
1622 [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]]
1623 [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]]
1624 [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]]
1625 [[1, 0, 2]] B? [[1, 0, 2]] B? [[1, 0, 2]]
1626 [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]]
1627 [[2, 1, 0]] B? [[2, 1, 0]] B? [[2, 1, 0]]
1628 [[0, 2, 1]] B? [[0, 2, 1]] B? [[0, 2, 1]]
1629 [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]]
1630 [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]] B? [[0, 2, 1], [1, 0, 2]]
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 [] BG [] BG []
1636 [] BG [] BG []
1637 [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]]
1638 [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]]
1639 [] BO [] BO []
1640 [] B_ [] B_ []
1641 [] BO [] BO []
1642 [] BO [] BO []
1643 [] BO [] BO []
1644 [] B_ [] B_ []
1645 [] BO [] BO []
1646 [] BO [] BO []
1647 [] BG [] BG []
1648 [] BG [] BG []
1649 [] BG [] BG []
1650 [[0, 2, 1]] B_ [[0, 2, 1]] B_ [[0, 2, 1]]
1651 [] BG [] BG []
1652 [[0, 2, 1]] B_ [[0, 2, 1]] B_ [[0, 2, 1]]
1653 [[0, 2, 1]] BG [[0, 2, 1]] BG [[0, 2, 1]]
1654 [[0, 2, 1]] BG [[0, 2, 1]] BG [[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 [] BO [] BO []
1660 [] B_ [] B_ []
1661 [] BO [] BO []
1662 [] BO [] BO []
1663 [] BG [] BG []
1664 [] BG [] BG []
1665 [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]]
1666 [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]]
1667 [] B_ [] B_ []
1668 [] BO [] BO []
1669 [] BO [] BO []
1670 [] BO [] BO []
1671 [] BG [] BG []
1672 [[2, 1, 0]] B_ [[2, 1, 0]] B_ [[2, 1, 0]]
1673 [] BG [] BG []
1674 [] BG [] BG []
1675 [[2, 1, 0]] B_ [[2, 1, 0]] B_ [[2, 1, 0]]
1676 [] BG [] BG []
1677 [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]]
1678 [[2, 1, 0]] BG [[2, 1, 0]] BG [[2, 1, 0]]
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 [] BW [] BW []
1684 [] Bg [] Bg []
1685 [] BW [] BW []
1686 [] BW [] BW []
1687 [] BW [] BW []
1688 [] Bg [] Bg []
1689 [] BW [] BW []
1690 [] BW [] BW []
1691 [] Bo [] Bo []
1692 [] Bo [] Bo []
1693 [[1, 0, 2]] Bo [[1, 0, 2]] Bo [[1, 0, 2]]
1694 [[1, 0, 2]] Bo [[1, 0, 2]] Bo [[1, 0, 2]]
1695 [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]]
1696 [] Bg [] Bg []
1697 [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]]
1698 [] Bg [] Bg []
1699 [] Bg [] Bg []
1700 [] Bg [] Bg []
1701 [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]]
1702 [[1, 0, 2]] BW [[1, 0, 2]] BW [[1, 0, 2]]
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 [] B_ [] B_ []
1708 [] BO [] BO []
1709 [] BO [] BO []
1710 [] BO [] BO []
1711 [] B_ [] B_ []
1712 [] BO [] BO []
1713 [] BO [] BO []
1714 [] BO [] BO []
1715 [] BG [] BG []
1716 [] BG [] BG []
1717 [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]]
1718 [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]]
1719 [[1, 0, 2]] B_ [[1, 0, 2]] B_ [[1, 0, 2]]
1720 [] BG [] BG []
1721 [[1, 0, 2]] B_ [[1, 0, 2]] B_ [[1, 0, 2]]
1722 [] BG [] BG []
1723 [] BG [] BG []
1724 [] BG [] BG []
1725 [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]]
1726 [[1, 0, 2]] BG [[1, 0, 2]] BG [[1, 0, 2]]
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 [] Bg [] Bg []
1732 [] BW [] BW []
1733 [] BW [] BW []
1734 [] BW [] BW []
1735 [] Bo [] Bo []
1736 [] Bo [] Bo []
1737 [[2, 1, 0]] Bo [[2, 1, 0]] Bo [[2, 1, 0]]
1738 [[2, 1, 0]] Bo [[2, 1, 0]] Bo [[2, 1, 0]]
1739 [] BW [] BW []
1740 [] Bg [] Bg []
1741 [] BW [] BW []
1742 [] BW [] BW []
1743 [] Bg [] Bg []
1744 [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]]
1745 [] Bg [] Bg []
1746 [] Bg [] Bg []
1747 [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]]
1748 [] Bg [] Bg []
1749 [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]]
1750 [[2, 1, 0]] BW [[2, 1, 0]] BW [[2, 1, 0]]
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 [] Bo [] Bo []
1756 [] Bo [] Bo []
1757 [[0, 2, 1]] Bo [[0, 2, 1]] Bo [[0, 2, 1]]
1758 [[0, 2, 1]] Bo [[0, 2, 1]] Bo [[0, 2, 1]]
1759 [] Bg [] Bg []
1760 [] BW [] BW []
1761 [] BW [] BW []
1762 [] BW [] BW []
1763 [] Bg [] Bg []
1764 [] BW [] BW []
1765 [] BW [] BW []
1766 [] BW [] BW []
1767 [] Bg [] Bg []
1768 [] Bg [] Bg []
1769 [] Bg [] Bg []
1770 [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]]
1771 [] Bg [] Bg []
1772 [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]]
1773 [[0, 2, 1]] BW [[0, 2, 1]] BW [[0, 2, 1]]
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 [] Bw [] Bw []
1780 [] Bw [] Bw []
1781 [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]]
1782 [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]]
1783 [] Bw [] Bw []
1784 [] Bw [] Bw []
1785 [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]]
1786 [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]]
1787 [] Bw [] Bw []
1788 [] Bw [] Bw []
1789 [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]]
1790 [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]]
1791 [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]]
1792 [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]]
1793 [[1, 0, 2]] Bw [[1, 0, 2]] Bw [[1, 0, 2]]
1794 [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]]
1795 [[2, 1, 0]] Bw [[2, 1, 0]] Bw [[2, 1, 0]]
1796 [[0, 2, 1]] Bw [[0, 2, 1]] Bw [[0, 2, 1]]
1797 [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]]
1798 [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]] Bw [[0, 2, 1], [1, 0, 2]]
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]]
1804 ::
1806 sage: C = graphs.CubeGraph(1)
1807 sage: gens = search_tree(C, [C.vertices()], lab=False)
1808 sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order()
1809 2
1810 sage: C = graphs.CubeGraph(2)
1811 sage: gens = search_tree(C, [C.vertices()], lab=False)
1812 sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order()
1813 8
1814 sage: C = graphs.CubeGraph(3)
1815 sage: gens = search_tree(C, [C.vertices()], lab=False)
1816 sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order()
1817 48
1818 sage: C = graphs.CubeGraph(4)
1819 sage: gens = search_tree(C, [C.vertices()], lab=False)
1820 sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order()
1821 384
1822 sage: C = graphs.CubeGraph(5)
1823 sage: gens = search_tree(C, [C.vertices()], lab=False)
1824 sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order()
1825 3840
1826 sage: C = graphs.CubeGraph(6)
1827 sage: gens = search_tree(C, [C.vertices()], lab=False)
1828 sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order()
1829 46080
1831 One can also turn off the indicator function (note- this will take
1832 longer)
1834 ::
1836 sage: D1 = DiGraph({0:[2],2:[0],1:[1]}, loops=True)
1837 sage: D2 = DiGraph({1:[2],2:[1],0:[0]}, loops=True)
1838 sage: a,b = search_tree(D1, [D1.vertices()], use_indicator_function=False)
1839 sage: c,d = search_tree(D2, [D2.vertices()], use_indicator_function=False)
1840 sage: b==d
1841 True
1843 Previously a bug, now the output is correct::
1845 sage: G = Graph('^????????????????????{??N??@w??FaGa?PCO@CP?AGa?_QO?Q@G?CcA??cc????Bo????{????F_')
1846 sage: perm = {3:15, 15:3}
1847 sage: H = G.relabel(perm, inplace=False)
1848 sage: G.canonical_label() == H.canonical_label()
1849 True
1851 Another former bug::
1853 sage: Graph('Fll^G').canonical_label()
1854 Graph on 7 vertices
1856 ::
1858 sage: g = Graph(21)
1859 sage: g.automorphism_group(return_group=False, order=True)
1860 51090942171709440000
1861 """
1862 cdef int i, j, m # local variables
1863 cdef int uif = 1 if use_indicator_function else 0
1864 cdef int _base = 1 if base else 0
1866 cdef OrbitPartition Theta
1867 cdef int index = 0 # see Theorem 2.33 in [1]
1868 size = Integer(1)
1870 cdef int L = 100 # memory limit for storing values from fix and mcr:
1871 # Phi and Omega store specific information about some
1872 # of the automorphisms we already know about, and they
1873 # are arrays of length L
1874 cdef int **Phi # stores the fixed point sets of each automorphism
1875 cdef int **Omega # stores the minimal elements of each cell of the
1876 # orbit partition
1877 cdef int l = -1 # current index for storing values in Phi and Omega-
1878 # we start at -1 so that when we increment first,
1879 # the first place we write to is 0.
1880 cdef int **W # for each k, W[k] is a list of the vertices to be searched
1881 # down from the current partition nest, at k
1882 # Phi and Omega are ultimately used to make the size of W
1883 # as small as possible
1885 cdef PartitionStack nu, zeta, rho
1886 cdef int k_rho # the number of partitions in rho
1887 cdef int h = -1 # longest common ancestor of zeta and nu:
1888 # zeta[h] == nu[h], zeta[h+1] != nu[h+1]
1889 cdef int hb # longest common ancestor of rho and nu:
1890 # rho[hb] == nu[hb], rho[hb+1] != nu[hb+1]
1891 cdef int hh = 1 # the height of the oldest ancestor of nu
1892 # satisfying Lemma 2.25 in [1]
1893 cdef int ht # smallest such that all descendants of zeta[ht]
1894 # are known to be equivalent
1896 cdef mpz_t *Lambda_mpz, *zf_mpz, *zb_mpz # for tracking indicator values
1897 # zf and zb are indicator vectors remembering Lambda[k] for zeta and rho,
1898 # respectively
1899 cdef int hzf # the max height for which Lambda and zf agree
1900 cdef int hzb = -1 # the max height for which Lambda and zb agree
1902 cdef CGraph M
1903 cdef int *gamma # for storing permutations
1904 cdef int *alpha # for storing pointers to cells of nu[k]:
1905 # allocated to be length 4*n + 1 for scratch (see functions
1906 # _sort_by_function and _refine)
1907 cdef int *v # list of vertices determining nu
1908 cdef int *e # 0 or 1, see states 12 and 17
1909 cdef int state # keeps track of place in algorithm
1910 cdef int _dig, tvh, n
1912 if isinstance(G, GenericGraph):
1913 n = G.order()
1914 elif isinstance(G, CGraph):
1915 M = G
1916 n = M.num_verts
1917 else:
1918 raise TypeError("G must be a Sage graph.")
1920 # trivial case
1921 if n == 0:
1922 if not (lab or dict_rep or certify or base or order):
1923 return [[]]
1924 output_tuple = [[[]]]
1925 if dict_rep:
1926 output_tuple.append({})
1927 if lab:
1928 output_tuple.append(G.copy())
1929 if certify:
1930 output_tuple.append({})
1931 if base:
1932 output_tuple.append([])
1933 if order:
1934 output_tuple.append(1)
1935 return tuple(output_tuple)
1937 # allocate int pointers
1938 W = <int **> sage_malloc( n * sizeof(int *) )
1939 Phi = <int **> sage_malloc( L * sizeof(int *) )
1940 Omega = <int **> sage_malloc( L * sizeof(int *) )
1942 # allocate GMP int pointers
1943 Lambda_mpz = <mpz_t *> sage_malloc( (n+2) * sizeof(mpz_t) )
1944 zf_mpz = <mpz_t *> sage_malloc( (n+2) * sizeof(mpz_t) )
1945 zb_mpz = <mpz_t *> sage_malloc( (n+2) * sizeof(mpz_t) )
1947 # check for memory errors
1948 if not (W and Phi and Omega and Lambda_mpz and zf_mpz and zb_mpz):
1949 sage_free(Lambda_mpz)
1950 sage_free(zf_mpz)
1951 sage_free(zb_mpz)
1952 sage_free(W)
1953 sage_free(Phi)
1954 sage_free(Omega)
1955 raise MemoryError("Error allocating memory.")
1957 # allocate int arrays
1958 gamma = <int *> sage_malloc( n * sizeof(int) )
1959 W[0] = <int *> sage_malloc( (n*n) * sizeof(int) )
1960 Phi[0] = <int *> sage_malloc( (L*n) * sizeof(int) )
1961 Omega[0] = <int *> sage_malloc( (L*n) * sizeof(int) )
1962 alpha = <int *> sage_malloc( (4*n + 1) * sizeof(int) )
1963 v = <int *> sage_malloc( n * sizeof(int) )
1964 e = <int *> sage_malloc( n * sizeof(int) )
1966 # check for memory errors
1967 if not (gamma and W[0] and Phi[0] and Omega[0] and alpha and v and e):
1968 sage_free(gamma)
1969 sage_free(W[0])
1970 sage_free(Phi[0])
1971 sage_free(Omega[0])
1972 sage_free(alpha)
1973 sage_free(v)
1974 sage_free(e)
1975 sage_free(Lambda_mpz)
1976 sage_free(zf_mpz)
1977 sage_free(zb_mpz)
1978 sage_free(W)
1979 sage_free(Phi)
1980 sage_free(Omega)
1981 raise MemoryError("Error allocating memory.")
1983 # setup double index arrays
1984 for i from 0 < i < n:
1985 W[i] = W[0] + n*i
1986 for i from 0 < i < L:
1987 Phi[i] = Phi[0] + n*i
1988 for i from 0 < i < L:
1989 Omega[i] = Omega[0] + n*i
1991 # allocate GMP ints
1992 for i from 0 <= i < n+2:
1993 mpz_init(Lambda_mpz[i])
1994 mpz_init_set_si(zf_mpz[i], -1) # correspond to default values of
1995 mpz_init_set_si(zb_mpz[i], -1) # "infinity"
1996 # Note that there is a potential memory leak here - if a particular
1997 # mpz fails to allocate, this is not checked for
1999 if isinstance(G, GenericGraph):
2000 # relabel vertices to the set {0,...,n-1}
2001 G = G.copy()
2002 ffrom = G.relabel(return_map=True)
2003 to = {}
2004 for vvv in ffrom.iterkeys():
2005 to[ffrom[vvv]] = vvv
2006 Pi2 = []
2007 for cell in Pi:
2008 Pi2.append([ffrom[c] for c in cell])
2009 Pi = Pi2
2010 if sparse:
2011 M = SparseGraph(n)
2012 else:
2013 M = DenseGraph(n)
2014 if isinstance(G, Graph):
2015 for i, j, la in G.edge_iterator():
2016 M.add_arc_unsafe(i,j)
2017 M.add_arc_unsafe(j,i)
2018 elif isinstance(G, DiGraph):
2019 for i, j, la in G.edge_iterator():
2020 M.add_arc_unsafe(i,j)
2022 # initialize W
2023 for i from 0 <= i < n:
2024 for j from 0 <= j < n:
2025 W[i][j] = 0
2027 # set up the rest of the variables
2028 nu = PartitionStack(Pi)
2029 Theta = OrbitPartition(n)
2030 output = []
2031 if dig: _dig = 1
2032 else: _dig = 0
2033 if certify:
2034 lab=True
2035 if _base:
2036 base_set = []
2038 if verbosity > 1:
2039 t = cputime()
2040 if verbosity > 2:
2041 rho = PartitionStack(n)
2042 zeta = PartitionStack(n)
2043 state = 1
2044 while state != -1:
2045 if verbosity > 0:
2046 print '-----'
2047 print 'k: ' + str(nu.k)
2048 print 'k_rho: ' + str(k_rho)
2049 print 'hh', hh
2050 print 'nu:'
2051 print nu
2052 if verbosity >= 2:
2053 t = cputime(t)
2054 print 'time:', t
2055 if verbosity >= 3:
2056 print 'zeta:'
2057 print [zeta.entries[iii] for iii in range(n)]
2058 print [zeta.levels[iii] for iii in range(n)]
2059 print 'rho'
2060 print [rho.entries[iii] for iii in range(n)]
2061 print [rho.levels[iii] for iii in range(n)]
2062 if verbosity >= 4:
2063 Thetarep = []
2064 for i from 0 <= i < n:
2065 j = Theta._find(i)
2066 didit = False
2067 for celll in Thetarep:
2068 if celll[0] == j:
2069 celll.append(i)
2070 didit = True
2071 if not didit:
2072 Thetarep.append([j])
2073 print 'Theta: ', str(Thetarep)
2074 if verbosity >= 5:
2075 if state == 1:
2076 verbose_first_time = True
2077 verbose_just_refined = False
2078 elif verbose_first_time:
2079 verbose_first_time = False
2080 # here we have just gone through step 1, and must now begin
2081 # to record information about the tree
2082 ST_vis = DiGraph()
2083 ST_vis_heights = {0:[nu.repr_at_k(0)]}
2084 ST_vis.add_vertex(nu.repr_at_k(0))
2085 ST_vis_current_vertex = nu.repr_at_k(0)
2086 ST_vis_current_level = 0
2087 #ST_vis.show(vertex_size=0)
2088 if state == 2:
2089 verbose_just_refined = True
2090 elif verbose_just_refined:
2091 verbose_just_refined = False
2092 # here we have gone through step 2, and must record the
2093 # refinement just made
2094 while ST_vis_current_level > nu.k-1:
2095 ST_vis_current_vertex = ST_vis.predecessors(ST_vis_current_vertex)[0]
2096 ST_vis_current_level -= 1
2097 if ST_vis_heights.has_key(nu.k):
2098 ST_vis_heights[nu.k].append(nu.repr_at_k(nu.k))
2099 else:
2100 ST_vis_heights[nu.k] = [nu.repr_at_k(nu.k)]
2101 ST_vis.add_edge(ST_vis_current_vertex, nu.repr_at_k(nu.k), '%d,%d'%(ST_vis_splitvert, ST_vis_invariant))
2102 ST_vis_current_vertex = nu.repr_at_k(nu.k)
2103 ST_vis_current_level += 1
2104 if state == 13 and nu.k == -1:
2105 ST_vis_new_heights = {}
2106 for ST_vis_k in ST_vis_heights:
2107 ST_vis_new_heights[-ST_vis_k] = ST_vis_heights[ST_vis_k]
2108 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()})
2109 print '-----'
2110 print 'state:', state
2113 if state == 1: # Entry point to algorithm
2114 # get alpha to point to cells of nu
2115 j = 1
2116 alpha[0] = 0
2117 for i from 0 < i < n:
2118 if nu.levels[i-1] == 0:
2119 alpha[j] = i
2120 j += 1
2121 alpha[j] = -1
2123 # "nu[0] := R(G, Pi, Pi)"
2124 nu._refine(alpha, n, M, _dig, uif)
2126 if not _dig:
2127 if nu._sat_225(n): hh = nu.k
2128 if nu._is_discrete(): state = 18; continue
2130 # store the first smallest nontrivial cell in W[k], and set v[k]
2131 # equal to its minimum element
2132 v[nu.k] = nu._first_smallest_nontrivial(W[nu.k], n)
2133 mpz_set_ui(Lambda_mpz[nu.k], 0)
2134 e[nu.k] = 0 # see state 12, and 17
2135 state = 2
2137 elif state == 2: # Move down the search tree one level by refining nu
2138 nu.k += 1
2140 # "nu[k] := nu[k-1] perp v[k-1]"
2141 nu._clear()
2142 alpha[0] = nu._split_vertex(v[nu.k-1])
2143 alpha[1] = -1
2144 i = nu._refine(alpha, n, M, _dig, uif)
2145 if verbosity >= 5:
2146 ST_vis_invariant = int(i)
2147 ST_vis_splitvert = int(v[nu.k-1])
2149 mpz_set_si(Lambda_mpz[nu.k], i)
2151 # only if this is the first time moving down the search tree:
2152 if h == -1: state = 5; continue
2154 # update hzf
2155 if hzf == nu.k-1 and mpz_cmp(Lambda_mpz[nu.k], zf_mpz[nu.k]) == 0: hzf = nu.k
2156 if not lab: state = 3; continue
2158 # "qzb := cmp(Lambda[k], zb[k])"
2159 if qzb == 0:
2160 if mpz_cmp_si(zb_mpz[nu.k], -1) == 0: # if "zb[k] == oo"
2161 qzb = -1
2162 else:
2163 qzb = mpz_cmp( Lambda_mpz[nu.k], zb_mpz[nu.k] )
2164 # update hzb
2165 if hzb == nu.k-1 and qzb == 0: hzb = nu.k
2167 # if Lambda[k] > zb[k], then zb[k] := Lambda[k]
2168 # (zb keeps track of the indicator invariants corresponding to
2169 # rho, the closest canonical leaf so far seen- if Lambda is
2170 # bigger, then rho must be about to change
2171 if qzb > 0: mpz_set(zb_mpz[nu.k], Lambda_mpz[nu.k])
2172 state = 3
2174 elif state == 3: # attempt to rule out automorphisms while moving down
2175 # the tree
2176 if hzf <= nu.k or (lab and qzb >= 0): # changed hzb to hzf, == to <=
2177 state = 4
2178 else: state = 6
2179 # if k > hzf, then we know that nu currently does not look like
2180 # zeta, the first terminal node encountered. Then if we are not
2181 # looking for a canonical label, there is no reason to continue.
2182 # However, if we are looking for one, and qzb < 0, i.e.
2183 # Lambda[k] < zb[k], then the indicator is not maximal, and we
2184 # can't reach a canonical leaf.
2186 elif state == 4: # at this point we have -not- ruled out the presence
2187 # of automorphisms
2188 if nu._is_discrete(): state = 7; continue
2190 # store the first smallest nontrivial cell in W[k], and set v[k]
2191 # equal to its minimum element
2192 v[nu.k] = nu._first_smallest_nontrivial(W[nu.k], n)
2194 if _dig or not nu._sat_225(n): hh = nu.k + 1
2195 e[nu.k] = 0 # see state 12, and 17
2196 state = 2 # continue down the tree
2198 elif state == 5: # alternative to 3: since we have not yet gotten
2199 # zeta, there are no automorphisms to rule out.
2200 # instead we record Lambda to zf and zb
2201 # (see state 3)
2202 if _base:
2203 base_set.append(v[nu.k-1])
2204 mpz_set(zf_mpz[nu.k], Lambda_mpz[nu.k])
2205 mpz_set(zb_mpz[nu.k], Lambda_mpz[nu.k])
2206 state = 4
2208 elif state == 6: # at this stage, there is no reason to continue
2209 # downward, and an automorphism has not been
2210 # discovered
2211 j = nu.k
2213 # return to the longest ancestor nu[i] of nu that could have a
2214 # descendant equivalent to zeta or could improve on rho.
2215 # All terminal nodes descending from nu[hh] are known to be
2216 # equivalent, so i < hh. Also, if i > hzb, none of the
2217 # descendants of nu[i] can improve rho, since the indicator is
2218 # off (Lambda(nu) < Lambda(rho)). If i >= ht, then no descendant
2219 # of nu[i] is equivalent to zeta (see [1, p67]).
2220 if ht-1 > hzb:
2221 if ht-1 < hh-1:
2222 nu.k = ht-1
2223 else:
2224 nu.k = hh-1
2225 else:
2226 if hzb < hh-1:
2227 nu.k = hzb
2228 else:
2229 nu.k = hh-1
2231 # TODO: investigate the following line
2232 if nu.k == -1: nu.k = 0 # not in BDM, broke at G = Graph({0:[], 1:[]}), Pi = [[0,1]], lab=False
2234 if lab:
2235 if hb > nu.k: # update hb since we are backtracking (not in [1])
2236 hb = nu.k # recall hb is the longest common ancestor of rho and nu
2238 if j == hh: state = 13; continue
2239 # recall hh: the height of the oldest ancestor of zeta for which
2240 # Lemma 2.25 is satsified, which implies that all terminal nodes
2241 # descended from there are equivalent (or simply k if 2.25 does
2242 # not apply). If we are looking at such a node, then the partition
2243 # at nu[hh] can be used for later pruning, so we store its fixed
2244 # set and a set of representatives of its cells
2245 if l < L-1: l += 1
2246 for i from 0 <= i < n:
2247 Omega[l][i] = 0 # changed Lambda to Omega
2248 Phi[l][i] = 0
2249 if nu._is_min_cell_rep(i, hh):
2250 Omega[l][i] = 1
2251 if nu._is_fixed(i, hh):
2252 Phi[l][i] = 1
2254 state = 12
2256 elif state == 7: # we have just arrived at a terminal node of the
2257 # search tree T(G, Pi)
2258 # if this is the first terminal node, go directly to 18, to
2259 # process zeta
2260 if h == -1: state = 18; continue
2262 # hzf is the extremal height of ancestors of both nu and zeta,
2263 # so if k < hzf, nu is not equivalent to zeta, i.e. there is no
2264 # automorphism to discover.
2265 # TODO: investigate why, in practice, the same does not seem to be
2266 # true for hzf < k... BDM had !=, not <, and this broke at
2267 # G = Graph({0:[],1:[],2:[]}), Pi = [[0,1,2]]
2268 if nu.k < hzf: state = 8; continue
2270 # get the permutation corresponding to this terminal node
2271 nu._get_permutation_from(zeta, gamma)
2273 if verbosity > 3:
2274 print 'checking for automorphism:'
2275 print [gamma[iii] for iii in range(n)]
2277 # if G^gamma == G, goto 10
2278 if _is_automorphism(M, n, gamma):
2279 state = 10
2280 else:
2281 state = 8
2283 elif state == 8: # we have just ruled out the presence of automorphism
2284 # and have not yet considered whether nu improves on
2285 # rho
2286 # if we are not searching for a canonical label, there is nothing
2287 # to do here
2288 if (not lab) or (qzb < 0):
2289 state = 6; continue
2291 # if Lambda[k] > zb[k] or nu is shorter than rho, then we have
2292 # found an improvement for rho
2293 if (qzb > 0) or (nu.k < k_rho):
2294 state = 9; continue
2296 # if G(nu) > G(rho) (returns 1), goto 9
2297 # if G(nu) < G(rho) (returns -1), goto 6
2298 # if G(nu) == G(rho) (returns 0), get the automorphism and goto 10
2299 m = nu._compare_with(M, n, rho)
2301 if m > 0:
2302 state = 9; continue
2303 if m < 0:
2304 state = 6; continue
2306 rho._get_permutation_from(nu, gamma)
2307 if verbosity > 3:
2308 print 'automorphism discovered:'
2309 print [gamma[iii] for iii in range(n)]
2310 state = 10
2312 elif state == 9: # entering this state, nu is a best-so-far guess at
2313 # the canonical label
2314 rho = PartitionStack(nu)
2315 k_rho = nu.k
2317 qzb = 0
2318 hb = nu.k
2319 hzb = nu.k
2321 # set zb[k+1] = Infinity
2322 mpz_set_si(zb_mpz[nu.k+1], -1)
2323 state = 6
2325 elif state == 10: # we have an automorphism to process
2326 # increment l
2327 if l < L - 1:
2328 l += 1
2330 for i from 0 <= i < n:
2331 if gamma[i] == i:
2332 Phi[l][i] = 1
2333 Omega[l][i] = 1
2334 else:
2335 Phi[l][i] = 0
2336 m = i
2337 j = gamma[i]
2338 while j != i:
2339 if j < m: m = j
2340 j = gamma[j]
2341 if m == i:
2342 Omega[l][i] = 1
2343 else:
2344 Omega[l][i] = 0
2345 m = Theta._incorporate_permutation(gamma, n)
2346 # if each orbit of gamma is part of an orbit in Theta, then the
2347 # automorphism is already in the span of those we have seen
2348 if not m:
2349 state = 11
2350 continue
2352 # record the automorphism
2353 output.append([ Integer(gamma[i]) for i from 0 <= i < n ])
2355 # The variable tvh represents the minimum element of W[k],
2356 # the last time we were at state 13 and backtracking up
2357 # zeta. If this is not still a minimal cell representative of Theta,
2358 # then we need to immediately backtrack to the place where it was
2359 # defined on a part of zeta, since the rest of the tree is now
2360 # equivalent. Otherwise, proceed to 11 and 12 before moving back to
2361 # the hub state.
2362 if Theta.elements[tvh] == -1:
2363 state = 11
2364 continue
2365 nu.k = h
2366 state = 13
2368 elif state == 11: # we have just discovered an automorphism,
2369 # but tvh is still a minimal representative for its
2370 # orbit in Theta. Therefore we cannot backtrack all
2371 # the way to where zeta meets nu. Instead we just use
2372 # indicator values to determine where to backtrack.
2373 if lab:
2374 nu.k = hb
2375 else:
2376 nu.k = h
2377 state = 12
2379 elif state == 12:
2380 # e keeps track of the whether W[k] has been thinned out by Omega
2381 # and Phi. It is set to 1 when you have just finished coming up the
2382 # search tree, and have intersected W[k] with Omega[i], for the
2383 # appropriate i < l, but since there may be an automorphism mapping
2384 # one element of W[k] to another still, we thin out W[k] again.
2385 # (see state 17) Coming from 11, this is an explicit automorphism.
2386 # Coming from 6, this is an implicit automorphism.
2387 if e[nu.k] == 1:
2388 for j from 0 <= j < n:
2389 if W[nu.k][j] and not Omega[l][j]:
2390 W[nu.k][j] = 0
2391 state = 13
2393 elif state == 13: # hub state
2394 if nu.k == -1:
2395 # the algorithm has finished
2396 state = -1; continue
2397 if nu.k > h:
2398 # if we are not at a node of zeta
2399 state = 17; continue
2400 if nu.k == h:
2401 # if we are at a node of zeta, then we have not yet backtracked
2402 # UP zeta, so skip the rest of 13
2403 state = 14; continue
2405 # thus, it must be that k < h, and this means we are done
2406 # searching underneath zeta[k+1], so now, k is the new longest
2407 # ancestor of nu and zeta:
2408 h = nu.k
2410 # set tvh to the minimum cell representative of W[k]
2411 # (see states 10 and 14)
2412 for i from 0 <= i < n:
2413 if W[nu.k][i]:
2414 tvh = i
2415 break
2416 state = 14
2418 elif state == 14: # iterate v[k] through W[k] until a minimum cell rep
2419 # of Theta is found:
2420 # this state gets hit only when we are looking for a
2421 # new split off of zeta
2422 # The variable tvh was set to be the minimum element of W[k]
2423 # the last time we were at state 13 and backtracking up
2424 # zeta. If this is in the same cell of Theta as v[k], increment
2425 # index (see Theorem 2.33 in [1])
2426 if Theta._find(v[nu.k]) == Theta._find(tvh):
2427 index += 1
2429 # find the next v[k] in W[k]
2430 i = v[nu.k] + 1
2431 while i < n and not W[nu.k][i]:
2432 i += 1
2433 if i < n:
2434 v[nu.k] = i
2435 else:
2436 # there is no new vertex to consider at this level
2437 v[nu.k] = -1
2438 state = 16
2439 continue
2441 # if the new v[k] is not a minimum cell representative of Theta,
2442 # then we already considered that rep., and that subtree was
2443 # isomorphic to the one corresponding to v[k]
2444 if Theta.elements[v[nu.k]] != -1: state = 14
2445 else:
2446 # otherwise, we do have a vertex to consider
2447 state = 15
2449 elif state == 15: # we have a new vertex, v[k], that we must split on
2450 # hh is smallest such that nu[hh] satisfies Lemma 2.25. If it is
2451 # larger than k+1, it must be modified, since we are changing that
2452 # part
2453 if nu.k + 1 < hh:
2454 hh = nu.k + 1
2455 # hzf is maximal such that indicators line up for nu and zeta
2456 if nu.k < hzf:
2457 hzf = nu.k
2458 if not lab or hzb < nu.k:
2459 # in either case there is no need to update hzb, which is the
2460 # length for which nu and rho have the same indicators
2461 state = 2; continue
2462 hzb = nu.k
2463 qzb = 0
2464 state = 2
2466 elif state == 16: # backtrack one level in the search tree, recording
2467 # information relevant to Theorem 2.33
2468 j = 0
2469 for i from 0 <= i < n:
2470 if W[nu.k][i]: j += 1
2471 if j == index and ht == nu.k+1: ht = nu.k
2472 size *= index
2473 index = 0
2474 nu.k -= 1
2476 if lab:
2477 if hb > nu.k: # update hb since we are backtracking (not in [1]):
2478 hb = nu.k # recall hb is the longest common ancestor of rho and nu
2480 state = 13
2482 elif state == 17: # you have just finished coming up the search tree,
2483 # and must now consider going back down.
2484 if e[nu.k] == 0:
2485 # intersect W[k] with each Omega[i] such that {v_0,...,v_(k-1)}
2486 # is contained in Phi[i]
2487 for i from 0 <= i <= l:
2488 # check if {v_0,...,v_(k-1)} is contained in Phi[i]
2489 # i.e. fixed pointwise by the automorphisms so far seen
2490 j = 0
2491 while j < nu.k and Phi[i][v[j]]:
2492 j += 1
2493 # if so, only check the minimal orbit representatives
2494 if j == nu.k:
2495 for j from 0 <= j < n:
2496 if W[nu.k][j] and not Omega[i][j]:
2497 W[nu.k][j] = 0
2498 e[nu.k] = 1 # see state 12
2500 # see if there is a relevant vertex to split on:
2501 i = v[nu.k] + 1
2502 while i < n and not W[nu.k][i]:
2503 i += 1
2504 if i < n:
2505 v[nu.k] = i
2506 state = 15
2507 continue
2508 else:
2509 v[nu.k] = -1
2511 # otherwise backtrack one level
2512 nu.k -= 1
2513 state = 13
2515 elif state == 18: # The first time we encounter a terminal node, we
2516 # come straight here to set up zeta. This is a one-
2517 # time state.
2518 # initialize counters for zeta:
2519 h = nu.k # zeta[h] == nu[h]
2520 ht = nu.k # nodes descended from zeta[ht] are all equivalent
2521 hzf = nu.k # max such that indicators for zeta and nu agree
2523 zeta = PartitionStack(nu)
2525 nu.k -= 1
2526 if not lab: state = 13; continue
2528 rho = PartitionStack(nu)
2530 # initialize counters for rho:
2531 k_rho = nu.k + 1 # number of partitions in rho
2532 hzb = nu.k # max such that indicators for rho and nu agree - BDM had k+1
2533 hb = nu.k # rho[hb] == nu[hb] - BDM had k+1
2535 qzb = 0 # Lambda[k] == zb[k], so...
2536 state = 13
2538 # free the GMP ints
2539 for i from 0 <= i < n+2:
2540 mpz_clear(Lambda_mpz[i])
2541 mpz_clear(zf_mpz[i])
2542 mpz_clear(zb_mpz[i])
2544 # free int arrays
2545 sage_free(gamma)
2546 sage_free(W[0])
2547 sage_free(Phi[0])
2548 sage_free(Omega[0])
2549 sage_free(alpha)
2550 sage_free(v)
2551 sage_free(e)
2553 # free GMP int pointers
2554 sage_free(Lambda_mpz)
2555 sage_free(zf_mpz)
2556 sage_free(zb_mpz)
2558 # free int pointers
2559 sage_free(W)
2560 sage_free(Phi)
2561 sage_free(Omega)
2563 # prepare output
2564 if not (lab or dict_rep or certify or base or order):
2565 return output
2566 output_tuple = [output]
2567 if dict_rep:
2568 if isinstance(G, GenericGraph):
2569 ddd = {}
2570 for vvv in ffrom.iterkeys(): # v is a C variable
2571 if ffrom[vvv] != 0:
2572 ddd[vvv] = ffrom[vvv]
2573 else:
2574 ddd[vvv] = n
2575 else: # G is a CGraph
2576 ddd = {}
2577 for i from 0 <= i < n:
2578 ddd[i] = i
2579 output_tuple.append(ddd)
2580 if lab:
2581 H = _term_pnest_graph(G, rho)
2582 output_tuple.append(H)
2583 if certify:
2584 dd = {}
2585 for i from 0 <= i < n:
2586 dd[to[rho.entries[i]]] = i
2587 output_tuple.append(dd)
2588 if base:
2589 output_tuple.append(base_set)
2590 if order:
2591 output_tuple.append(size)
2592 return tuple(output_tuple)
2594 # Benchmarking functions
2596 def all_labeled_graphs(n):
2597 """
2598 Returns all labeled graphs on n vertices 0,1,...,n-1. Used in
2599 classifying isomorphism types (naive approach), and more
2600 importantly in benchmarking the search algorithm.
2602 EXAMPLE::
2604 sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree
2605 sage: from sage.graphs.graph_isom import all_labeled_graphs
2606 sage: Glist = {}
2607 sage: Giso = {}
2608 sage: for n in range(1,5):
2609 ... Glist[n] = all_labeled_graphs(n)
2610 ... Giso[n] = []
2611 ... for g in Glist[n]:
2612 ... a, b = search_tree(g, [range(n)])
2613 ... inn = False
2614 ... for gi in Giso[n]:
2615 ... if b == gi:
2616 ... inn = True
2617 ... if not inn:
2618 ... Giso[n].append(b)
2619 sage: for n in Giso:
2620 ... print n, len(Giso[n])
2621 1 1
2622 2 2
2623 3 4
2624 4 11
2625 sage: n = 5
2626 sage: Glist[n] = all_labeled_graphs(n)
2627 sage: Giso[n] = []
2628 sage: for g in Glist[5]:
2629 ... a, b = search_tree(g, [range(n)])
2630 ... inn = False
2631 ... for gi in Giso[n]:
2632 ... if b == gi:
2633 ... inn = True
2634 ... if not inn:
2635 ... Giso[n].append(b)
2636 sage: print n, len(Giso[n]) # long time
2637 5 34
2638 sage: graphs_list.show_graphs(Giso[4])
2639 """
2640 TE = []
2641 for i in range(n):
2642 for j in range(i):
2643 TE.append((i, j))
2644 m = len(TE)
2645 Glist= []
2646 for i in range(2**m):
2647 G = Graph(n)
2648 b = Integer(i).binary()
2649 b = '0'*(m-len(b)) + b
2650 for i in range(m):
2651 if int(b[i]):
2652 G.add_edge(TE[i])
2653 Glist.append(G)
2654 return Glist
2656 def kpow(listy, k):
2657 """
2658 Returns the subset of the power set of listy consisting of subsets
2659 of size k. Used in all_ordered_partitions.
2661 EXAMPLE::
2663 sage: from sage.graphs.graph_isom import kpow
2664 sage: kpow(['a', 1, {}], 2)
2665 [[1, 'a'], [{}, 'a'], ['a', 1], [{}, 1], ['a', {}], [1, {}]]
2666 """
2667 list = []
2668 if k > 1:
2669 for LL in kpow(listy, k-1):
2670 for a in listy:
2671 if not a in LL:
2672 list.append([a] + LL)
2673 if k == 1:
2674 for i in listy:
2675 list.append([i])
2676 return list
2678 def all_ordered_partitions(listy):
2679 """
2680 Returns all ordered partitions of the set 0,1,...,n-1. Used in
2681 benchmarking the search algorithm.
2683 EXAMPLE::
2685 sage: from sage.graphs.graph_isom import all_ordered_partitions
2686 sage: all_ordered_partitions(['a', 1, {}])
2687 [[['a'], [1], [{}]],
2688 [['a'], [{}], [1]],
2689 [['a'], [{}, 1]],
2690 [['a'], [1, {}]],
2691 [[1], ['a'], [{}]],
2692 [[1], [{}], ['a']],
2693 [[1], [{}, 'a']],
2694 [[1], ['a', {}]],
2695 [[{}], ['a'], [1]],
2696 [[{}], [1], ['a']],
2697 [[{}], [1, 'a']],
2698 [[{}], ['a', 1]],
2699 [[1, 'a'], [{}]],
2700 [[{}, 'a'], [1]],
2701 [['a', 1], [{}]],
2702 [[{}, 1], ['a']],
2703 [['a', {}], [1]],
2704 [[1, {}], ['a']],
2705 [[{}, 1, 'a']],
2706 [[1, {}, 'a']],
2707 [[{}, 'a', 1]],
2708 [['a', {}, 1]],
2709 [[1, 'a', {}]],
2710 [['a', 1, {}]]]
2711 """
2712 LL = []
2713 for i in range(1,len(listy)+1):
2714 for cell in kpow(listy, i):
2715 list_remainder = [x for x in listy if x not in cell]
2716 remainder_partitions = all_ordered_partitions(list_remainder)
2717 for remainder in remainder_partitions:
2718 LL.append( [cell] + remainder )
2719 if len(listy) == 0:
2720 return [[]]
2721 else:
2722 return LL
2724 def all_labeled_digraphs_with_loops(n):
2725 """
2726 Returns all labeled digraphs on n vertices 0,1,...,n-1. Used in
2727 classifying isomorphism types (naive approach), and more
2728 importantly in benchmarking the search algorithm.
2730 EXAMPLE::
2732 sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree
2733 sage: from sage.graphs.graph_isom import all_labeled_digraphs_with_loops
2734 sage: Glist = {}
2735 sage: Giso = {}
2736 sage: for n in range(1,4):
2737 ... Glist[n] = all_labeled_digraphs_with_loops(n)
2738 ... Giso[n] = []
2739 ... for g in Glist[n]:
2740 ... a, b = search_tree(g, [range(n)], dig=True)
2741 ... inn = False
2742 ... for gi in Giso[n]:
2743 ... if b == gi:
2744 ... inn = True
2745 ... if not inn:
2746 ... Giso[n].append(b)
2747 sage: for n in Giso:
2748 ... print n, len(Giso[n])
2749 1 2
2750 2 10
2751 3 104
2752 """
2753 TE = []
2754 for i in range(n):
2755 for j in range(n):
2756 TE.append((i, j))
2757 m = len(TE)
2758 Glist= []
2759 for i in range(2**m):
2760 G = DiGraph(n, loops=True)
2761 b = Integer(i).binary()
2762 b = '0'*(m-len(b)) + b
2763 for j in range(m):
2764 if int(b[j]):
2765 G.add_edge(TE[j])
2766 Glist.append(G)
2767 return Glist
2769 def all_labeled_digraphs(n):
2770 """
2771 EXAMPLES::
2773 sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import search_tree
2774 sage: from sage.graphs.graph_isom import all_labeled_digraphs
2775 sage: Glist = {}
2776 sage: Giso = {}
2777 sage: for n in range(1,4):
2778 ... Glist[n] = all_labeled_digraphs(n)
2779 ... Giso[n] = []
2780 ... for g in Glist[n]:
2781 ... a, b = search_tree(g, [range(n)], dig=True)
2782 ... inn = False
2783 ... for gi in Giso[n]:
2784 ... if b == gi:
2785 ... inn = True
2786 ... if not inn:
2787 ... Giso[n].append(b)
2788 sage: for n in Giso:
2789 ... print n, len(Giso[n])
2790 1 1
2791 2 3
2792 3 16
2793 """
2795 ## sage: n = 4 # long time (4 minutes)
2796 ## sage: Glist[n] = all_labeled_digraphs(n) # long time
2797 ## sage: Giso[n] = [] # long time
2798 ## sage: for g in Glist[n]: # long time
2799 ## ... a, b = search_tree(g, [range(n)], dig=True)
2800 ## ... inn = False
2801 ## ... for gi in Giso[n]:
2802 ## ... if enum(b) == enum(gi):
2803 ## ... inn = True
2804 ## ... if not inn:
2805 ## ... Giso[n].append(b)
2806 ## sage: print n, len(Giso[n]) # long time
2807 ## 4 218
2808 TE = []
2809 for i in range(n):
2810 for j in range(n):
2811 if i != j: TE.append((i, j))
2812 m = len(TE)
2813 Glist= []
2814 for i in range(2**m):
2815 G = DiGraph(n, loops=True)
2816 b = Integer(i).binary()
2817 b = '0'*(m-len(b)) + b
2818 for j in range(m):
2819 if int(b[j]):
2820 G.add_edge(TE[j])
2821 Glist.append(G)
2822 return Glist
2824 def perm_group_elt(lperm):
2825 """
2826 Given a list permutation of the set 0, 1, ..., n-1, returns the
2827 corresponding PermutationGroupElement where we take 0 = n.
2829 EXAMPLE::
2831 sage: from sage.graphs.graph_isom import perm_group_elt
2832 sage: perm_group_elt([0,2,1])
2833 (1,2)
2834 sage: perm_group_elt([1,2,0])
2835 (1,2,3)
2836 """
2837 from sage.groups.perm_gps.permgroup_named import SymmetricGroup
2838 n = len(lperm)
2839 S = SymmetricGroup(n)
2840 Part = orbit_partition(lperm, list_perm=True)
2841 gens = []
2842 for z in Part:
2843 if len(z) > 1:
2844 if 0 in z:
2845 zed = z.index(0)
2846 generator = z[:zed] + [n] + z[zed+1:]
2847 gens.append(tuple(generator))
2848 else:
2849 gens.append(tuple(z))
2850 E = S(gens)
2851 return E
2853 def orbit_partition(gamma, list_perm=False):
2854 r"""
2855 Assuming that G is a graph on vertices 0,1,...,n-1, and gamma is an
2856 element of SymmetricGroup(n), returns the partition of the vertex
2857 set determined by the orbits of gamma, considered as action on the
2858 set 1,2,...,n where we take 0 = n. In other words, returns the
2859 partition determined by a cyclic representation of gamma.
2861 INPUT:
2864 - ``list_perm`` - if True, assumes
2865 ``gamma`` is a list representing the map
2866 `i \mapsto ``gamma``[i]`.
2869 EXAMPLES::
2871 sage: from sage.graphs.graph_isom import orbit_partition
2872 sage: G = graphs.PetersenGraph()
2873 sage: S = SymmetricGroup(10)
2874 sage: gamma = S('(10,1,2,3,4)(5,6,7)(8,9)')
2875 sage: orbit_partition(gamma)
2876 [[1, 2, 3, 4, 0], [5, 6, 7], [8, 9]]
2877 sage: gamma = S('(10,5)(1,6)(2,7)(3,8)(4,9)')
2878 sage: orbit_partition(gamma)
2879 [[1, 6], [2, 7], [3, 8], [4, 9], [5, 0]]
2880 """
2881 if list_perm:
2882 n = len(gamma)
2883 seen = [1] + [0]*(n-1)
2884 i = 0
2885 p = 0
2886 partition = [[0]]
2887 while sum(seen) < n:
2888 if gamma[i] != partition[p][0]:
2889 partition[p].append(gamma[i])
2890 i = gamma[i]
2891 seen[i] = 1
2892 else:
2893 i = min([j for j in range(n) if seen[j] == 0])
2894 partition.append([i])
2895 p += 1
2896 seen[i] = 1
2897 return partition
2898 else:
2899 n = len(gamma.list())
2900 l = []
2901 for i in range(1,n+1):
2902 orb = gamma.orbit(i)
2903 if orb not in l: l.append(orb)
2904 for i in l:
2905 for j in range(len(i)):
2906 if i[j] == n:
2907 i[j] = 0
2908 return l
2910 def verify_partition_refinement(G, initial_partition, terminal_partition):
2911 """
2912 Verify that the refinement is correct.
2914 EXAMPLE::
2916 sage: from sage.graphs.graph_isom import PartitionStack
2917 sage: from sage.graphs.base.sparse_graph import SparseGraph
2918 sage: G = SparseGraph(10)
2919 sage: for i,j,_ in graphs.PetersenGraph().edge_iterator():
2920 ... G.add_arc(i,j)
2921 ... G.add_arc(j,i)
2922 sage: P = PartitionStack(10)
2923 sage: P.set_k(1)
2924 sage: P.split_vertex(0)
2925 sage: P.refine(G, [0], 10, 0, 1)
2926 sage: P
2927 (0,2,3,6,7,8,9,1,4,5)
2928 (0|2,3,6,7,8,9|1,4,5)
2929 sage: P.set_k(2)
2930 sage: P.split_vertex(1)
2932 Note that this line implicitly tests the function
2933 verify_partition_refinement::
2935 sage: P.refine(G, [7], 10, 0, 1, test=True)
2936 sage: P
2937 (0,3,7,8,9,2,6,1,4,5)
2938 (0|3,7,8,9,2,6|1,4,5)
2939 (0|3,7,8,9|2,6|1|4,5)
2940 """
2941 if not G.is_equitable(terminal_partition):
2942 raise RuntimeError("Resulting partition is not equitable!!!!!!!!!\n"+\
2943 str(initial_partition) + "\n" + \
2944 str(terminal_partition) + "\n" + \
2945 str(G.am()))