Ticket #6112: trac_6112.patch

File trac_6112.patch, 115.0 KB (added by ddrake, 10 years ago)
  • deleted file sage/graphs/graph_isom.pyx

    # HG changeset patch
    # User Dan Drake <drake@kaist.edu>
    # Date 1245118851 -32400
    # Node ID 24f862891b4b242cc8222c5e048a247fe271c297
    # Parent  d22a81385eb6030f5904e7c0604722b0ac305b81
    [mq]: trac_6112
    
    diff --git a/sage/graphs/graph_isom.pyx b/sage/graphs/graph_isom.pyx
    deleted file mode 100644
    + -  
    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