Ticket #3676: trac3676-refactor_graph_isom.patch

sage/graphs/graph_generators.py

@@ -2987 +2987 @@
         import networkx
         return graph.DiGraph(networkx.gnc_graph(n, seed))
 
+    def RandomDirectedGNP(self, n, p):
+        r"""
+        Returns a random digraph on $n$ nodes.  Each edge is inserted
+        independently with probability $p$.
+
+        REFERENCES:
+            [1] P. Erdos and A. Renyi, On Random Graphs, Publ. Math. 6, 290 (1959).
+            [2] E. N. Gilbert, Random Graphs, Ann. Math. Stat., 30, 1141 (1959).
+        
+        PLOTTING:
+        When plotting, this graph will use the default spring-layout
+        algorithm, unless a position dictionary is specified.
+                
+        EXAMPLE:
+            sage: digraphs.RandomDirectedGNP(10, .2).num_verts()
+            10
+        
+        """
+        from random import random
+        D = graph.DiGraph(n)
+        for i in xrange(n):
+            for j in xrange(i):
+                if random() < p:
+                    D.add_edge(i,j)
+            for j in xrange(i+1,n):
+                if random() < p:
+                    D.add_edge(i,j)
+        return D
+
     def RandomDirectedGNR(self, n, p, seed=None):
         """
         Returns a random GNR (growing network with redirection) digraph with n

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/__init__.py

@@ +1 @@
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

@@ +1 @@
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include '../../../ext/cdefs.pxi'
+include '../../../ext/stdsage.pxi'
+
+from sage.graphs.base.c_graph cimport CGraph
+from sage.graphs.base.sparse_graph cimport SparseGraph
+from sage.graphs.base.dense_graph cimport DenseGraph
+from sage.rings.integer cimport Integer
+from sage.groups.perm_gps.partn_ref.tree_traversal cimport PartitionStack, PS_is_discrete, PS_move_min_to_front, PS_num_cells, traverse_tree, traverse_return
+
+cdef class GraphStruct:
+    cdef CGraph G
+    cdef bint directed
+    cdef bint use_indicator
+    cdef int *scratch # length 3n+1
+
+cdef int refine_by_degree(PartitionStack *, object, int *, int)
+cdef int compare_graphs(int *, int *, object)
+cdef bint all_children_are_equivalent(PartitionStack *, object)
+cdef int degree(PartitionStack *, CGraph, int, int, bint)
+cdef int sort_by_function(PartitionStack *, int, int *)
+
+
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

@@ +1 @@
+"""
+Module for graph-theoretic partition backtrack functions.
+
+DOCTEST:
+    sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
+
+"""
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+def search_tree(G_in, partition, lab=True, dig=False, dict_rep=False, certify=False,
+                    verbosity=0, use_indicator_function=True, sparse=False,
+                    base=False, order=False):
+    """
+    Compute canonical labels and automorphism groups of graphs.
+    
+    INPUT:
+    G_in -- a Sage graph
+    partition -- a list of lists representing a partition of the vertices
+    lab -- if True, return the canonical label in addition to the automorphism
+        group.
+    dig -- if True, does not use Lemma 2.25 in [1], and the algorithm is valid
+        for digraphs and graphs with loops.
+    dict_rep -- if True, explain which vertices are which elements of the set
+        {1,2,...,n} in the representation of the automorphism group.
+    certify -- if True, return the relabeling from G to its canonical
+        label.
+    verbosity -- currently ignored
+    use_indicator_function -- option to turn off indicator function
+        (False -> slower)
+    sparse -- whether to use sparse or dense representation of the graph
+        (ignored if G is already a CGraph - see sage.graphs.base)
+    base -- whether to return the first sequence of split vertices (used in
+        computing the order of the group)
+    order -- whether to return the order of the automorphism group
+
+    OUTPUT:
+    Depends on the options. If more than one thing is returned, they are in a
+    tuple in the following order:
+        list of generators in list-permutation format -- always
+        dict -- if dict_rep
+        graph -- if lab
+        dict -- if certify
+        list -- if base
+        integer -- if order
+
+    DOCTEST:
+        sage: st = sage.groups.perm_gps.partn_ref.refinement_graphs.search_tree
+        sage: from sage.graphs.base.dense_graph import DenseGraph
+        sage: from sage.graphs.base.sparse_graph import SparseGraph
+
+    Graphs on zero vertices:
+        sage: G = Graph()
+        sage: st(G, [[]], order=True)
+        ([], Graph on 0 vertices, 1)
+
+    Graphs on one vertex:
+        sage: G = Graph(1)
+        sage: st(G, [[0]], order=True)
+        ([], Graph on 1 vertex, 1)
+
+    Graphs on two vertices:
+        sage: G = Graph(2)
+        sage: st(G, [[0,1]], order=True)
+        ([[1, 0]], Graph on 2 vertices, 2)
+        sage: st(G, [[0],[1]], order=True)
+        ([], Graph on 2 vertices, 1)
+        sage: G.add_edge(0,1)
+        sage: st(G, [[0,1]], order=True)
+        ([[1, 0]], Graph on 2 vertices, 2)
+        sage: st(G, [[0],[1]], order=True)
+        ([], Graph on 2 vertices, 1)
+
+    Graphs on three vertices:
+        sage: G = Graph(3)
+        sage: st(G, [[0,1,2]], order=True)
+        ([[0, 2, 1], [1, 0, 2]], Graph on 3 vertices, 6)
+        sage: st(G, [[0],[1,2]], order=True)
+        ([[0, 2, 1]], Graph on 3 vertices, 2)
+        sage: st(G, [[0],[1],[2]], order=True)
+        ([], Graph on 3 vertices, 1)
+        sage: G.add_edge(0,1)
+        sage: st(G, [range(3)], order=True)
+        ([[1, 0, 2]], Graph on 3 vertices, 2)
+        sage: st(G, [[0],[1,2]], order=True)
+        ([], Graph on 3 vertices, 1)
+        sage: st(G, [[0,1],[2]], order=True)
+        ([[1, 0, 2]], Graph on 3 vertices, 2)
+
+    The Dodecahedron has automorphism group of size 120:
+        sage: G = graphs.DodecahedralGraph()
+        sage: Pi = [range(20)]
+        sage: st(G, Pi, order=True)[2]
+        120
+
+    The three-cube has automorphism group of size 48:
+        sage: G = graphs.CubeGraph(3)
+        sage: G.relabel()
+        sage: Pi = [G.vertices()]
+        sage: st(G, Pi, order=True)[2]
+        48
+
+    We obtain the same output using different types of Sage graphs:
+        sage: G = graphs.DodecahedralGraph()
+        sage: GD = DenseGraph(20)
+        sage: GS = SparseGraph(20)
+        sage: for i,j,_ in G.edge_iterator():
+        ...    GD.add_arc(i,j); GD.add_arc(j,i)
+        ...    GS.add_arc(i,j); GS.add_arc(j,i)
+        sage: Pi=[range(20)]
+        sage: a,b = st(G, Pi)
+        sage: asp,bsp = st(GS, Pi)
+        sage: ade,bde = st(GD, Pi)
+        sage: bsg = Graph(implementation='networkx')
+        sage: bdg = Graph(implementation='networkx')
+        sage: for i in range(20):
+        ...    for j in range(20):
+        ...        if bsp.has_arc(i,j):
+        ...            bsg.add_edge(i,j)
+        ...        if bde.has_arc(i,j):
+        ...            bdg.add_edge(i,j)
+        sage: print a, b.graph6_string()
+        [[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
+        sage: a == asp
+        True
+        sage: a == ade
+        True
+        sage: b == bsg
+        True
+        sage: b == bdg
+        True
+
+    Cubes!
+        sage: C = graphs.CubeGraph(1)
+        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
+        2
+        sage: C = graphs.CubeGraph(2)
+        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
+        8
+        sage: C = graphs.CubeGraph(3)
+        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
+        48
+        sage: C = graphs.CubeGraph(4)
+        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
+        384
+        sage: C = graphs.CubeGraph(5)
+        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
+        3840
+        sage: C = graphs.CubeGraph(6)
+        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
+        46080
+
+    One can also turn off the indicator function (note- this will take longer)
+        sage: D1 = DiGraph({0:[2],2:[0],1:[1]}, loops=True)
+        sage: D2 = DiGraph({1:[2],2:[1],0:[0]}, loops=True)
+        sage: a,b = st(D1, [D1.vertices()], dig=True, use_indicator_function=False)
+        sage: c,d = st(D2, [D2.vertices()], dig=True, use_indicator_function=False)
+        sage: b==d
+        True
+
+    This example is due to Chris Godsil:
+        sage: HS = graphs.HoffmanSingletonGraph()
+        sage: clqs = (HS.complement()).cliques()
+        sage: alqs = [Set(c) for c in clqs if len(c) == 15]
+        sage: Y = Graph([alqs, lambda s,t: len(s.intersection(t))==0], implementation='networkx')
+        sage: Y0,Y1 = Y.connected_components_subgraphs()
+        sage: st(Y0, [Y0.vertices()])[1] == st(Y1, [Y1.vertices()])[1]
+        True
+        sage: st(Y0, [Y0.vertices()])[1] == st(HS, [HS.vertices()])[1]
+        True
+        sage: st(HS, [HS.vertices()])[1] == st(Y1, [Y1.vertices()])[1]
+        True
+
+    Certain border cases need to be tested as well:
+        sage: G = Graph('Fll^G')
+        sage: a,b,c = st(G, [range(G.num_verts())], order=True); b
+        Graph on 7 vertices
+        sage: c
+        48
+        sage: G = Graph(21)
+        sage: st(G, [range(G.num_verts())], order=True)[2] == factorial(21)
+        True
+
+        sage: G = Graph('^????????????????????{??N??@w??FaGa?PCO@CP?AGa?_QO?Q@G?CcA??cc????Bo????{????F_')
+        sage: perm = {3:15, 15:3}
+        sage: H = G.relabel(perm, inplace=False)
+        sage: st(G, [range(G.num_verts())])[1] == st(H, [range(H.num_verts())])[1]
+        True
+
+    """
+    cdef CGraph G
+    cdef int i, j, n
+    cdef Integer I
+    cdef traverse_return *output
+    cdef int **part
+    from sage.graphs.graph import GenericGraph, Graph, DiGraph
+    if isinstance(G_in, GenericGraph):
+        n = G_in.num_verts()
+        G_in = G_in.copy()
+        if G_in.vertices() != range(n):
+            to = G_in.relabel(return_map=True)
+            frm = {}
+            for v in to.iterkeys():
+                frm[to[v]] = v
+            partition = [[to[v] for v in cell] for cell in partition]
+        else:
+            to = range(n)
+            frm = to
+        if sparse:
+            G = SparseGraph(n)
+        else:
+            G = DenseGraph(n)
+        if G_in.is_directed():
+            for i from 0 <= i < n:
+                for _,j,_ in G_in.outgoing_edge_iterator(i):
+                    G.add_arc(i,j)
+        else:
+            for i from 0 <= i < n:
+                for _,j,_ in G_in.edge_iterator(i):
+                    if j <= i:
+                        G.add_arc(i,j)
+                        G.add_arc(j,i)
+    elif isinstance(G_in, CGraph):
+        G = <CGraph> G_in
+        n = G.num_verts
+        to = range(n)
+        frm = to
+    else:
+        raise TypeError("G must be a Sage graph.")
+    
+    part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
+    if part is NULL:
+        raise MemoryError
+    for i from 0 <= i < len(partition):
+        part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
+        if part[i] is NULL:
+            raise MemoryError
+        for j from 0 <= j < len(partition[i]):
+            part[i][j] = partition[i][j]
+        part[i][len(partition[i])] = -1
+    part[len(partition)] = NULL
+
+    cdef GraphStruct GS = GraphStruct()
+    GS.G = G
+    GS.directed = 1 if dig else 0
+    GS.use_indicator = 1 if use_indicator_function else 0
+    GS.scratch = <int *> sage_malloc( (3*G.num_verts + 1) * sizeof(int) )
+    if GS.scratch is NULL:
+        raise MemoryError
+
+    output = traverse_tree(GS, part, G.num_verts, &all_children_are_equivalent, &refine_by_degree, &compare_graphs, lab, base, order)
+    sage_free( GS.scratch )
+    # prepare output
+    list_of_gens = []
+    for i from 0 <= i < output.num_gens:
+        list_of_gens.append([output.generators[j+i*G.num_verts] for j from 0 <= j < G.num_verts])
+    return_tuple = [list_of_gens]
+    if dict_rep:
+        ddd = {}
+        for v in frm.iterkeys():
+            if frm[v] != 0:
+                ddd[v] = frm[v]
+            else:
+                ddd[v] = n
+        return_tuple.append(ddd)
+    if lab:
+        if isinstance(G_in, GenericGraph):
+            G_C = G_in.copy()
+            G_C.relabel([output.relabeling[i] for i from 0 <= i < n])
+        else:
+            if isinstance(G, SparseGraph):
+                G_C = SparseGraph(n)
+            else:
+                G_C = DenseGraph(n)
+            for i from 0 <= i < n:
+                for j in G.out_neighbors(i):
+                    G_C.add_arc(output.relabeling[i],output.relabeling[j])
+        return_tuple.append(G_C)
+    if certify:
+        dd = {}
+        for i from 0 <= i < G.num_verts:
+            dd[frm[i]] = output.relabeling[i]
+        return_tuple.append(dd)
+    if base:
+        return_tuple.append([output.base[i] for i from 0 <= i < output.base_size])
+    if order:
+        I = Integer()
+        mpz_set(I.value, output.order)
+        return_tuple.append(I)
+    for i from 0 <= i < len(partition):
+        sage_free(part[i])
+    sage_free(part)
+    mpz_clear(output.order)
+    sage_free(output.generators)
+    if base:
+        sage_free(output.base)
+    if lab:
+        sage_free(output.relabeling)
+    sage_free(output)
+    if len(return_tuple) == 1:
+        return return_tuple[0]
+    else:
+        return tuple(return_tuple)
+
+cdef int refine_by_degree(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
+    """
+    Refines the input partition by checking degrees of vertices to the given
+    cells.
+    
+    INPUT:
+    PS -- a partition stack, whose finest partition is the partition to be
+        refined.
+    S -- a graph struct object, which contains scratch space, the graph in
+        question, and some flags.
+    cells_to_refine_by -- a list of pointers to cells to check degrees against
+        in refining the other cells (updated in place)
+    ctrb_len -- how many cells in the above
+    
+    OUTPUT:
+    An invariant, such that if the function is called with a permutation of the
+    situation (i.e. there is a gamma taking one graph to the other, which also
+    takes one partition to the other), the same number is returned.
+    """
+    cdef GraphStruct GS = <GraphStruct> S
+    cdef CGraph G = GS.G
+    cdef int current_cell_against = 0
+    cdef int current_cell, i, r
+    cdef int first_largest_subcell
+    cdef int invariant = 1
+    cdef int *degrees = GS.scratch # length 3n+1
+    cdef bint necessary_to_split_cell
+    cdef int against_index
+    while not PS_is_discrete(PS) and current_cell_against < ctrb_len:
+        invariant += 1
+        current_cell = 0
+        while current_cell < PS.degree:
+            invariant += 50
+            i = current_cell
+            necessary_to_split_cell = 0
+            while 1:
+                degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 0)
+                if degrees[i-current_cell] != degrees[0]:
+                    necessary_to_split_cell = 1
+                i += 1
+                if PS.levels[i-1] <= PS.depth:
+                    break
+            # now, i points to the next cell (before refinement)
+            if necessary_to_split_cell:
+                invariant += 10
+                first_largest_subcell = sort_by_function(PS, current_cell, degrees)
+                invariant += first_largest_subcell
+                against_index = current_cell_against
+                while against_index < ctrb_len:
+                    if cells_to_refine_by[against_index] == current_cell:
+                        cells_to_refine_by[against_index] = first_largest_subcell
+                        break
+                    against_index += 1
+                against_index = ctrb_len
+                r = current_cell
+                while 1:
+                    if r == current_cell or PS.levels[r-1] == PS.depth:
+                        if r != first_largest_subcell:
+                            cells_to_refine_by[against_index] = r
+                            against_index += 1
+                            ctrb_len += 1
+                    r += 1
+                    if r >= i:
+                        break
+                invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 0)
+                invariant += (i - current_cell)
+            current_cell = i
+        if GS.directed:
+            # if we are looking at a digraph, also compute
+            # the reverse degrees and sort by them
+            current_cell = 0
+            while current_cell < PS.degree: # current_cell is still a valid cell
+                invariant += 20
+                i = current_cell
+                necessary_to_split_cell = 0
+                while 1:
+                    degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 1)
+                    if degrees[i-current_cell] != degrees[0]:
+                        necessary_to_split_cell = 1
+                    i += 1
+                    if PS.levels[i-1] <= PS.depth:
+                        break
+                # now, i points to the next cell (before refinement)
+                if necessary_to_split_cell:
+                    invariant += 7
+                    first_largest_subcell = sort_by_function(PS, current_cell, degrees)
+                    invariant += first_largest_subcell
+                    against_index = current_cell_against
+                    while against_index < ctrb_len:
+                        if cells_to_refine_by[against_index] == current_cell:
+                            cells_to_refine_by[against_index] = first_largest_subcell
+                            break
+                        against_index += 1
+                    against_index = ctrb_len
+                    r = current_cell
+                    while 1:
+                        if r == current_cell or PS.levels[r-1] == PS.depth:
+                            if r != first_largest_subcell:
+                                cells_to_refine_by[against_index] = r
+                                against_index += 1
+                                ctrb_len += 1
+                        r += 1
+                        if r >= i:
+                            break
+                    invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 1)
+                    invariant += (i - current_cell)
+                current_cell = i
+        current_cell_against += 1
+    if GS.use_indicator:
+        return invariant
+    else:
+        return 0
+    
+cdef int compare_graphs(int *gamma_1, int *gamma_2, object S):
+    """
+    Compare gamma_1(G) and gamma_2(G), returning 0 iff gamma_1(G) == gamma_2(G).
+    
+    INPUT:
+    gamma_1, gamma_2 -- list permutations (inverse)
+    S -- a graph struct object
+    
+    """
+    cdef int i, j
+    cdef GraphStruct GS = <GraphStruct> S
+    cdef CGraph G = GS.G
+    for i from 0 <= i < G.num_verts:
+        for j from 0 <= j < G.num_verts:
+            if G.has_arc_unsafe(gamma_1[i], gamma_1[j]):
+                if not G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
+                    return 1
+            elif G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
+                return -1
+    return 0
+
+cdef bint all_children_are_equivalent(PartitionStack *PS, object S):
+    """
+    Returns True if any refinement of the current partition results in the same
+    structure. Converse does not hold in general.
+    
+    INPUT:
+    PS -- the partition stack to be checked
+    S -- a graph struct object
+    """
+    cdef GraphStruct GS = <GraphStruct> S
+    cdef CGraph G = GS.G
+    cdef int i, n = PS.degree
+    cdef bint in_cell = 0
+    cdef int nontrivial_cells = 0
+    cdef int total_cells = PS_num_cells(PS)
+    if n <= total_cells + 4:
+        return 1
+    for i from 0 <= i < n-1:
+        if PS.levels[i] <= PS.depth:
+            if in_cell:
+                nontrivial_cells += 1
+            in_cell = 0
+        else:
+            in_cell = 1
+    if in_cell:
+        nontrivial_cells += 1
+    if n == total_cells + nontrivial_cells:
+        return 1
+    if n == total_cells + nontrivial_cells + 1:
+        return 1
+    return 0
+
+cdef int degree(PartitionStack *PS, CGraph G, int entry, int cell_index, bint reverse):
+    """
+    Returns the number of edges from the vertex corresponding to entry to
+    vertices in the cell corresponding to cell_index.
+    
+    INPUT:
+    PS -- the partition stack to be checked
+    S -- a graph struct object
+    entry -- the position of the vertex in question in the entries of PS
+    cell_index -- the starting position of the cell in question in the entries
+        of PS
+    reverse -- whether to check for arcs in the other direction
+    """
+    cdef int num_arcs = 0
+    entry = PS.entries[entry]
+    if not reverse:
+        while 1:
+            if G.has_arc_unsafe(PS.entries[cell_index], entry):
+                num_arcs += 1
+            if PS.levels[cell_index] > PS.depth:
+                cell_index += 1
+            else:
+                break
+    else:
+        while 1:
+            if G.has_arc_unsafe(entry, PS.entries[cell_index]):
+                num_arcs += 1
+            if PS.levels[cell_index] > PS.depth:
+                cell_index += 1
+            else:
+                break
+    return num_arcs
+
+cdef int sort_by_function(PartitionStack *PS, int start, int *degrees):
+    """
+    A simple counting sort, given the degrees of vertices to a certain cell.
+    
+    INPUT:
+    PS -- the partition stack to be checked
+    start -- beginning index of the cell to be sorted
+    degrees -- the values to be sorted by
+
+    """
+    cdef int n = PS.degree
+    cdef int i, j, m = 2*n, max, max_location
+    cdef int *counts = degrees + n, *output = degrees + 2*n + 1
+    for i from 0 <= i <= n:
+        counts[i] = 0
+    i = 0
+    while PS.levels[i+start] > PS.depth:
+        counts[degrees[i]] += 1
+        i += 1
+    counts[degrees[i]] += 1
+    # i+start is the right endpoint of the cell now
+    max = counts[0]
+    max_location = 0
+    for j from 0 < j <= n:
+        if counts[j] > max:
+            max = counts[j]
+            max_location = j
+        counts[j] += counts[j - 1]
+    for j from i >= j >= 0:
+        counts[degrees[j]] -= 1
+        output[counts[degrees[j]]] = PS.entries[start+j]
+    max_location = counts[max_location]+start
+    for j from 0 <= j <= i:
+        PS.entries[start+j] = output[j]
+    j = 1
+    while j <= n and counts[j] <= i:
+        if counts[j] > 0:
+            PS.levels[start + counts[j] - 1] = PS.depth
+        PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1)
+        j += 1
+    return max_location
+    
+    
+def all_labeled_graphs(n):
+    """
+    Returns all labeled graphs on n vertices {0,1,...,n-1}. Used in
+    classifying isomorphism types (naive approach), and more importantly
+    in benchmarking the search algorithm.
+    
+    EXAMPLE:
+        sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import all_labeled_graphs
+        sage: st = sage.groups.perm_gps.partn_ref.refinement_graphs.search_tree
+        sage: Glist = {}
+        sage: Giso  = {}
+        sage: for n in [1..5]:
+        ...    Glist[n] = all_labeled_graphs(n)
+        ...    Giso[n] = []
+        ...    for g in Glist[n]:
+        ...        a, b = st(g, [range(n)])
+        ...        inn = False
+        ...        for gi in Giso[n]:
+        ...            if b == gi:
+        ...                inn = True
+        ...        if not inn:
+        ...            Giso[n].append(b)
+        sage: for n in Giso:
+        ...    print n, len(Giso[n])
+        1 1
+        2 2
+        3 4
+        4 11
+        5 34
+
+    """
+    from sage.graphs.graph import Graph
+    TE = []
+    for i in range(n):
+        for j in range(i):
+            TE.append((i, j))
+    m = len(TE)
+    Glist= []
+    for i in range(2**m):
+        G = Graph(n)
+        b = Integer(i).binary()
+        b = '0'*(m-len(b)) + b
+        for i in range(m):
+            if int(b[i]):
+                G.add_edge(TE[i])
+        Glist.append(G)
+    return Glist
+
+
+def random_tests(t=10.0, n_max=60, perms_per_graph=10):
+    """
+    Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma
+    and random graphs G.
+    
+    DOCTEST:
+        sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
+        sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests()
+        All passed: ... random tests on ... graphs.
+        sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests(180.0, 200, 30) # long time
+        All passed: ... random tests on ... graphs.        
+
+    """
+    from sage.misc.misc import walltime
+    from sage.misc.prandom import random, randint
+    from sage.graphs.graph_generators import GraphGenerators, DiGraphGenerators
+    from sage.combinat.permutation import Permutations
+    cdef int i, j, num_tests = 0, num_graphs = 0, passed = 1
+    GG = GraphGenerators()
+    DGG = DiGraphGenerators()
+    t_0 = walltime()
+    while walltime(t_0) < t:
+        p = random()
+        n = randint(1, n_max)
+        S = Permutations(n)
+
+        G = GG.RandomGNP(n, p)
+        H = G.copy()
+        for i from 0 <= i < perms_per_graph:
+            G1 = search_tree(G, [G.vertices()])[1]
+            perm = list(S.random_element())
+            perm = [perm[j]-1 for j from 0 <= j < n]
+            G.relabel(perm)
+            G2 = search_tree(G, [G.vertices()])[1]
+            if G1 != G2:
+                print "FAILURE: graph6-"
+                print H.graph6_string()
+                print perm
+                passed = 0
+
+        D = DGG.RandomDirectedGNP(n, p)
+        D.loops(True)
+        for i from 0 <= i < n:
+            if random() < p:
+                D.add_edge(i,i)
+        E = D.copy()
+        for i from 0 <= i < perms_per_graph:
+            D1 = search_tree(D, [D.vertices()], dig=True)[1]
+            perm = list(S.random_element())
+            perm = [perm[j]-1 for j from 0 <= j < n]
+            D.relabel(perm)
+            D2 = search_tree(D, [D.vertices()], dig=True)[1]
+            if D1 != D2:
+                print "FAILURE: dig6-"
+                print E.dig6_string()
+                print perm
+                passed = 0
+        num_tests += 2*perms_per_graph
+        num_graphs += 2
+    if passed:
+        print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs)
+
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/tree_traversal.pxd

@@ +1 @@
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include '../../../ext/cdefs.pxi'
+include '../../../ext/stdsage.pxi'
+from sage.rings.integer cimport Integer
+
+cdef struct OrbitPartition:
+    # Disjoint set data structure for representing the orbits of the generators
+    # so far found. Also keeps track of the minimum elements of the cells and
+    # their sizes.
+    int degree
+    int *parent
+    int *rank
+    int *mcr # minimum cell representatives - only valid at the root of a cell
+    int *size # also only valid at the root of a cell
+
+cdef inline OrbitPartition *OP_new(int)
+cdef inline int OP_dealloc(OrbitPartition *)
+cdef inline int OP_find(OrbitPartition *, int)
+cdef inline int OP_join(OrbitPartition *, int, int)
+cdef inline int OP_merge_list_perm(OrbitPartition *, int *)
+
+cdef struct PartitionStack:
+    # Representation of a node of the search tree. A sequence of partitions of
+    # length depth + 1, each of which is finer than the last. Partition k is
+    # represented as PS.entries in order, broken up immediately after each
+    # entry of levels which is at most k.
+    int *entries
+    int *levels
+    int depth
+    int degree
+
+cdef inline PartitionStack *PS_new(int, bint)
+cdef inline PartitionStack *PS_copy(PartitionStack *)
+cdef inline int PS_copy_from_to(PartitionStack *, PartitionStack *)
+cdef inline PartitionStack *PS_from_list(object)
+cdef inline int PS_dealloc(PartitionStack *)
+cdef inline int PS_print(PartitionStack *)
+cdef inline int PS_print_partition(PartitionStack *, int)
+cdef inline bint PS_is_discrete(PartitionStack *)
+cdef inline int PS_num_cells(PartitionStack *)
+cdef inline int PS_move_min_to_front(PartitionStack *, int, int)
+cdef inline bint PS_is_mcr(PartitionStack *, int)
+cdef inline bint PS_is_fixed(PartitionStack *, int)
+cdef inline int PS_clear(PartitionStack *)
+cdef inline int PS_split_point(PartitionStack *, int)
+#cdef inline int PS_first_smallest(PartitionStack *, bitset_t)
+cdef inline int PS_get_perm_from(PartitionStack *, PartitionStack *, int *)
+
+cdef inline int split_point_and_refine(PartitionStack *, int, object,
+    int (*)(PartitionStack *, object, int *, int), int *)
+
+cdef struct traverse_return:
+    int *generators
+    int num_gens
+    int *relabeling
+    int *base
+    int base_size
+    mpz_t order
+
+cdef traverse_return *traverse_tree( object, int **, int,
+    bint (*)(PartitionStack *, object),
+    int (*)(PartitionStack *, object, int *, int),
+    int (*)(int *, int *, object), bint, bint, bint)
+
+
+

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/tree_traversal.pyx

@@ +1 @@
+"""
+Module for traversing the search tree associated with partition backtrack
+problems.
+
+DOCTEST:
+    sage: import sage.groups.perm_gps.partn_ref.tree_traversal
+
+"""
+
+#*****************************************************************************
+#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
+#
+# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
+#                         http://www.gnu.org/licenses/
+#*****************************************************************************
+
+include '../../../misc/bitset.pxi'
+
+cdef inline int min(int a, int b):
+    if a < b: return a
+    else: return b
+
+cdef inline int max(int a, int b):
+    if a > b: return a
+    else: return b
+
+cdef inline int cmp(int a, int b):
+    if a < b: return -1
+    elif a == b: return 0
+    else: return 1
+
+# OrbitPartitions
+
+cdef inline OrbitPartition *OP_new(int n):
+    """
+    Allocate and return a pointer to a new OrbitPartition of degree n. Returns a
+    null pointer in the case of an allocation failure.
+    """
+    cdef int i
+    cdef OrbitPartition *OP = <OrbitPartition *> \
+                                sage_malloc(sizeof(OrbitPartition))
+    if OP is NULL:
+        return OP
+    OP.degree = n
+    OP.parent = <int *> sage_malloc( n * sizeof(int) )
+    OP.rank = <int *> sage_malloc( n * sizeof(int) )
+    OP.mcr = <int *> sage_malloc( n * sizeof(int) )
+    OP.size = <int *> sage_malloc( n * sizeof(int) )
+    if OP.parent is NULL or OP.rank is NULL or \
+       OP.mcr is NULL or OP.size is NULL:
+        if OP.parent is not NULL: sage_free(OP.parent)
+        if OP.rank is not NULL: sage_free(OP.rank)
+        if OP.mcr is not NULL: sage_free(OP.mcr)
+        if OP.size is not NULL: sage_free(OP.size)
+        return NULL
+    for i from 0 <= i < n:
+        OP.parent[i] = i
+        OP.rank[i] = 0
+        OP.mcr[i] = i
+        OP.size[i] = 1
+    return OP
+
+cdef inline int OP_dealloc(OrbitPartition *OP):
+    sage_free(OP.parent)
+    sage_free(OP.rank)
+    sage_free(OP.mcr)
+    sage_free(OP.size)
+    sage_free(OP)
+    return 0
+
+cdef inline int OP_find(OrbitPartition *OP, int n):
+    """
+    Report the representative ("root") of the cell which contains n.
+    """
+    if OP.parent[n] == n:
+        return n
+    else:
+        OP.parent[n] = OP_find(OP, OP.parent[n])
+        return OP.parent[n]
+
+cdef inline int OP_join(OrbitPartition *OP, int m, int n):
+    """
+    Join the cells containing m and n, if they are different.
+    """
+    cdef int m_root = OP_find(OP, m)
+    cdef int n_root = OP_find(OP, n)
+    if OP.rank[m_root] > OP.rank[n_root]:
+        OP.parent[n_root] = m_root
+        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
+        OP.size[m_root] += OP.size[n_root]
+    elif OP.rank[m_root] < OP.rank[n_root]:
+        OP.parent[m_root] = n_root
+        OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root])
+        OP.size[n_root] += OP.size[m_root]
+    elif m_root != n_root:
+        OP.parent[n_root] = m_root
+        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
+        OP.size[m_root] += OP.size[n_root]
+        OP.rank[m_root] += 1
+
+cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma):
+    """
+    Joins the cells of OP which intersect the same orbit of gamma.
+    
+    INPUT:
+        gamma - an integer array representing i -> gamma[i].
+    
+    OUTPUT:
+        1 - something changed
+        0 - orbits of gamma all contained in cells of OP
+    """
+    cdef int i, i_root, gamma_i_root, changed = 0
+    for i from 0 <= i < OP.degree:
+        if gamma[i] == i: continue
+        i_root = OP_find(OP, i)
+        gamma_i_root = OP_find(OP, gamma[i])
+        if i_root != gamma_i_root:
+            changed = 1
+            OP_join(OP, i_root, gamma_i_root)
+    return changed
+
+def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]):
+    """
+    Demonstration and testing.
+    
+    DOCTEST:
+        sage: from sage.groups.perm_gps.partn_ref.tree_traversal \
+        ... import OP_represent
+        sage: OP_represent()
+        Allocating OrbitPartition...
+        Allocation passed.
+        Checking that each element reports itself as its root.
+        Each element reports itself as its root.
+        Merging:
+        Merged 0 and 1.
+        Merged 2 and 3.
+        Merged 3 and 4.
+        Done merging.
+        Finding:
+        0 -> 0, root: size=2, mcr=0, rank=1
+        1 -> 0
+        2 -> 2, root: size=3, mcr=2, rank=1
+        3 -> 2
+        4 -> 2
+        5 -> 5, root: size=1, mcr=5, rank=0
+        6 -> 6, root: size=1, mcr=6, rank=0
+        7 -> 7, root: size=1, mcr=7, rank=0
+        8 -> 8, root: size=1, mcr=8, rank=0
+        Allocating array to test merge_perm.
+        Allocation passed.
+        Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8]
+        Done merging.
+        Finding:
+        0 -> 0, root: size=5, mcr=0, rank=2
+        1 -> 0
+        2 -> 0
+        3 -> 0
+        4 -> 0
+        5 -> 5, root: size=3, mcr=5, rank=1
+        6 -> 5
+        7 -> 5
+        8 -> 8, root: size=1, mcr=8, rank=0
+        Deallocating OrbitPartition.
+        Done.
+
+    """
+    cdef int i
+    print "Allocating OrbitPartition..."
+    cdef OrbitPartition *OP = OP_new(n)
+    if OP is NULL:
+        print "Allocation failed!"
+        return
+    print "Allocation passed."
+    print "Checking that each element reports itself as its root."
+    good = True
+    for i from 0 <= i < n:
+        if not OP_find(OP, i) == i:
+            print "Failed at i = %d!"%i
+            good = False
+    if good: print "Each element reports itself as its root."
+    print "Merging:"
+    for i,j in merges:
+        OP_join(OP, i, j)
+        print "Merged %d and %d."%(i,j)
+    print "Done merging."
+    print "Finding:"
+    for i from 0 <= i < n:
+        j = OP_find(OP, i)
+        s = "%d -> %d"%(i, j)
+        if i == j:
+            s += ", root: size=%d, mcr=%d, rank=%d"%\
+                   (OP.size[i], OP.mcr[i], OP.rank[i])
+        print s
+    print "Allocating array to test merge_perm."
+    cdef int *gamma = <int *> sage_malloc( n * sizeof(int) )
+    if gamma is NULL:
+        print "Allocation failed!"
+        OP_dealloc(OP)
+        return
+    print "Allocation passed."
+    for i from 0 <= i < n:
+        gamma[i] = perm[i]
+    print "Merging permutation: %s"%perm
+    OP_merge_list_perm(OP, gamma)
+    print "Done merging."
+    print "Finding:"
+    for i from 0 <= i < n:
+        j = OP_find(OP, i)
+        s = "%d -> %d"%(i, j)
+        if i == j:
+            s += ", root: size=%d, mcr=%d, rank=%d"%\
+                   (OP.size[i], OP.mcr[i], OP.rank[i])
+        print s    
+    print "Deallocating OrbitPartition."
+    sage_free(gamma)
+    OP_dealloc(OP)
+    print "Done."
+
+# PartitionStacks
+
+cdef inline PartitionStack *PS_new(int n, bint unit_partition):
+    """
+    Allocate and return a pointer to a new PartitionStack of degree n. Returns a
+    null pointer in the case of an allocation failure.
+    """
+    cdef int i
+    cdef PartitionStack *PS = <PartitionStack *> \
+                                sage_malloc(sizeof(PartitionStack))
+    if PS is NULL:
+        return PS
+    PS.entries = <int *> sage_malloc( n * sizeof(int) )
+    PS.levels = <int *> sage_malloc( n * sizeof(int) )
+    if PS.entries is NULL or PS.levels is NULL:
+        if PS.entries is not NULL: sage_free(PS.entries)
+        if PS.levels is not NULL: sage_free(PS.levels)
+        return NULL
+    PS.depth = 0
+    PS.degree = n
+    if unit_partition:
+        for i from 0 <= i < n-1:
+            PS.entries[i] = i
+            PS.levels[i] = n
+        PS.entries[n-1] = n-1
+        PS.levels[n-1] = -1
+    return PS
+
+cdef inline PartitionStack *PS_copy(PartitionStack *PS):
+    """
+    Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
+    pointer in the case of an allocation failure.
+    """
+    cdef int i
+    cdef PartitionStack *PS2 = <PartitionStack *> \
+                                sage_malloc(sizeof(PartitionStack))
+    if PS2 is NULL:
+        return PS2
+    PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
+    PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
+    if PS2.entries is NULL or PS2.levels is NULL:
+        if PS2.entries is not NULL: sage_free(PS2.entries)
+        if PS2.levels is not NULL: sage_free(PS2.levels)
+        return NULL
+    PS_copy_from_to(PS, PS2)
+    return PS2
+
+cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2):
+    """
+    Copy all data from PS to PS2.
+    """
+    PS2.depth = PS.depth
+    PS2.degree = PS.degree
+    for i from 0 <= i < PS.degree:
+        PS2.entries[i] = PS.entries[i]
+        PS2.levels[i] = PS.levels[i]    
+
+cdef inline PartitionStack *PS_from_list(object L):
+    """
+    Allocate and return a pointer to a PartitionStack representing L. Returns a
+    null pointer in the case of an allocation failure.
+    """
+    cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0
+    for cell from 0 <= cell < num_cells:
+        n += len(L[cell])
+    cdef PartitionStack *PS = PS_new(n, False)
+    cell = 0
+    if PS is NULL:
+        return PS
+    while 1:
+        cur_len = len(L[cell])
+        for i from 0 <= i < cur_len:
+            PS.entries[cur_start + i] = L[cell][i]
+            PS.levels[cur_start + i] = n
+        PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
+        cur_start += cur_len
+        cell += 1
+        if cell == num_cells:
+            PS.levels[cur_start-1] = -1
+            break
+        PS.levels[cur_start-1] = 0
+    return PS
+
+cdef inline int PS_dealloc(PartitionStack *PS):
+    sage_free(PS.entries)
+    sage_free(PS.levels)
+    sage_free(PS)
+    return 0
+
+cdef inline int PS_print(PartitionStack *PS):
+    """
+    Print a visual representation of PS.
+    """
+    cdef int i
+    for i from 0 <= i < PS.depth + 1:
+        PS_print_partition(PS, i)
+
+cdef inline int PS_print_partition(PartitionStack *PS, int k):
+    """
+    Print the partition at depth k.
+    """
+    s = '('
+    for i from 0 <= i < PS.degree:
+        s += str(PS.entries[i])
+        if PS.levels[i] <= k:
+            s += '|'
+        else:
+            s += ','
+    s = s[:-1] + ')'
+    print s
+
+cdef inline bint PS_is_discrete(PartitionStack *PS):
+    """
+    Returns whether the deepest partition consists only of singleton cells.
+    """
+    cdef int i
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] > PS.depth:
+            return 0
+    return 1
+
+cdef inline int PS_num_cells(PartitionStack *PS):
+    """
+    Returns the number of cells.
+    """
+    cdef int i, ncells = 0
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] <= PS.depth:
+            ncells += 1
+    return ncells
+
+cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
+    """
+    Makes sure that the first element of the segment of entries i with
+    start <= i <= end is minimal.
+    """
+    cdef int i, min_loc = start, minimum = PS.entries[start]
+    for i from start < i <= end:
+        if PS.entries[i] < minimum:
+            min_loc = i
+            minimum = PS.entries[i]
+    if min_loc != start:
+        PS.entries[min_loc] = PS.entries[start]
+        PS.entries[start] = minimum
+
+cdef inline bint PS_is_mcr(PartitionStack *PS, int m):
+    """
+    Returns whether PS.elements[m] (not m!) is the smallest element of its cell.
+    """
+    return m == 0 or PS.levels[m-1] <= PS.depth
+
+cdef inline bint PS_is_fixed(PartitionStack *PS, int m):
+    """
+    Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming
+    PS_is_mcr(PS, m) is already true.
+    """
+    return PS.levels[m] <= PS.depth
+
+cdef inline int PS_clear(PartitionStack *PS):
+    """
+    Sets the current partition to the first shallower one, i.e. forgets about
+    boundaries between cells that are new to the current level.
+    """
+    cdef int i, cur_start = 0
+    for i from 0 <= i < PS.degree:
+        if PS.levels[i] >= PS.depth:
+            PS.levels[i] += 1
+        if PS.levels[i] < PS.depth:
+            PS_move_min_to_front(PS, cur_start, i)
+            cur_start = i+1
+
+cdef inline int PS_split_point(PartitionStack *PS, int v):
+    """
+    Detaches the point v from the cell it is in, putting the singleton cell
+    of just v in front. Returns the position where v is now located.
+    """
+    cdef int i = 0, index_of_v
+    while PS.entries[i] != v:
+        i += 1
+    index_of_v = i
+    while PS.levels[i] > PS.depth:
+        i += 1
+    if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \
+       and PS.levels[index_of_v] > PS.depth:
+        # if v is first (make sure v is not already alone)
+        PS_move_min_to_front(PS, index_of_v+1, i)
+        PS.levels[index_of_v] = PS.depth
+        return index_of_v
+    else:
+        # If v is not at front, i.e. v is not minimal in its cell,
+        # then move_min_to_front is not necessary since v will swap
+        # with the first before making its own cell, leaving it at
+        # the front of the other.
+        i = index_of_v
+        while i != 0 and PS.levels[i-1] > PS.depth:
+            i -= 1
+        PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v
+        PS.entries[i+1] = PS.entries[i] # move the first (min) to second
+        PS.entries[i] = v # place v first
+        PS.levels[i] = PS.depth
+        return i
+
+cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
+    """
+    Find the first occurrence of the smallest cell of size greater than one,
+    store its entries to b, and return its minimum element.
+    """
+    cdef int i = 0, j = 0, location = 0, n = PS.degree
+    bitset_zero(b)
+    while 1:
+        if PS.levels[i] <= PS.depth:
+            if i != j and n > i - j + 1:
+                n = i - j + 1
+                location = j
+            j = i + 1
+        if PS.levels[i] == -1: break
+        i += 1
+    # location now points to the beginning of the first, smallest,
+    # nontrivial cell
+    i = location
+    while 1:
+        bitset_flip(b, PS.entries[i])
+        if PS.levels[i] <= PS.depth: break
+        i += 1
+    return PS.entries[location]
+
+cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
+    """
+    Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
+    denotes the entry of the ith cell of the discrete partition PS.
+    """
+    cdef int i = 0
+    for i from 0 <= i < PS1.degree:
+        gamma[PS2.entries[i]] = PS1.entries[i]
+
+def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]):
+    """
+    Demonstration and testing.
+    
+    DOCTEST:
+        sage: from sage.groups.perm_gps.partn_ref.tree_traversal \
+        ... import PS_represent
+        sage: PS_represent()
+        Allocating PartitionStack...
+        Allocation passed:
+        (0,1,2,3,4,5,6,7,8,9)
+        Checking that entries are in order and correct level.
+        Everything seems in order, deallocating.
+        Deallocated.
+        Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]].
+        PartitionStack's data:
+        entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2]
+        levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1]
+        depth = 0, degree = 10
+        (6|3,4,8,7|1,9,5|0,2)
+        Checking PS_is_discrete:
+        False
+        Checking PS_num_cells:
+        4
+        Checking PS_is_mcr, min cell reps are:
+        [6, 3, 1, 0]
+        Checking PS_is_fixed, fixed elements are:
+        [6]
+        Copying PartitionStack:
+        (6|3,4,8,7|1,9,5|0,2)
+        Checking for consistency.
+        Everything is consistent.
+        Clearing copy:
+        (0,3,4,8,7,1,9,5,6,2)
+        Splitting point 6 from original:
+        0
+        (6|3,4,8,7|1,9,5|0,2)
+        Splitting point 1 from original:
+        5
+        (6|3,4,8,7|1|5,9|0,2)
+        Splitting point 8 from original:
+        1
+        (6|8|3,4,7|1|5,9|0,2)
+        Splitting point 2 from original:
+        8
+        (6|8|3,4,7|1|5,9|2|0)
+        Getting permutation from PS2->PS:
+        [6, 1, 0, 8, 3, 9, 2, 7, 4, 5]
+        Finding first smallest:
+        Minimal element is 5, bitset is:
+        0000010001
+        Deallocating PartitionStacks.
+        Done.
+
+    """
+    cdef int i, n = sum([len(cell) for cell in partition]), *gamma
+    cdef bitset_t b
+    print "Allocating PartitionStack..."
+    cdef PartitionStack *PS = PS_new(n, True), *PS2
+    if PS is NULL:
+        print "Allocation failed!"
+        return
+    print "Allocation passed:"
+    PS_print(PS)
+    print "Checking that entries are in order and correct level."
+    good = True
+    for i from 0 <= i < n-1:
+        if not (PS.entries[i] == i and PS.levels[i] == n):
+            print "Failed at i = %d!"%i
+            print PS.entries[i], PS.levels[i], i, n
+            good = False
+    if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1):
+        print "Failed at i = %d!"%(n-1)
+        good = False
+    if not PS.degree == n or not PS.depth == 0:
+        print "Incorrect degree or depth!"
+        good = False
+    if good: print "Everything seems in order, deallocating."
+    PS_dealloc(PS)
+    print "Deallocated."
+    print "Creating PartitionStack from partition %s."%partition
+    PS = PS_from_list(partition)
+    print "PartitionStack's data:"
+    print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n]
+    print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n]
+    print "depth = %d, degree = %d"%(PS.depth,PS.degree)
+    PS_print(PS)
+    print "Checking PS_is_discrete:"
+    print "True" if PS_is_discrete(PS) else "False"
+    print "Checking PS_num_cells:"
+    print PS_num_cells(PS)
+    print "Checking PS_is_mcr, min cell reps are:"
+    L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)]
+    print L
+    print "Checking PS_is_fixed, fixed elements are:"
+    print [PS.entries[l] for l in L if PS_is_fixed(PS, l)]
+    print "Copying PartitionStack:"
+    PS2 = PS_copy(PS)
+    PS_print(PS2)
+    print "Checking for consistency."
+    good = True
+    for i from 0 <= i < n:
+        if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]:
+            print "Failed at i = %d!"%i
+            good = False
+    if PS.degree != PS2.degree or PS.depth != PS2.depth:
+        print "Failure with degree or depth!"
+        good = False
+    if good:
+        print "Everything is consistent."
+    print "Clearing copy:"
+    PS_clear(PS2)
+    PS_print(PS2)
+    for s in splits:
+        print "Splitting point %d from original:"%s
+        print PS_split_point(PS, s)
+        PS_print(PS)
+    print "Getting permutation from PS2->PS:"
+    gamma = <int *> sage_malloc(n * sizeof(int))
+    PS_get_perm_from(PS, PS2, gamma)
+    print [gamma[i] for i from 0 <= i < n]
+    sage_free(gamma)
+    print "Finding first smallest:"
+    bitset_init(b, n)
+    i = PS_first_smallest(PS, b)
+    print "Minimal element is %d, bitset is:"%i
+    print bitset_string(b)
+    bitset_clear(b)
+    print "Deallocating PartitionStacks."
+    PS_dealloc(PS)
+    PS_dealloc(PS2)
+    print "Done."
+
+# Functions
+
+cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
+    int (*refine_and_return_invariant)\
+         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
+    int *cells_to_refine_by):
+    """
+    Make the partition stack one longer by copying the last partition in the
+    stack, split off a given point, and refine. Return the invariant given by
+    the refinement function.
+    
+    INPUT:
+    PS -- the partition stack to refine
+    v -- the point to split
+    S -- the structure
+    refine_and_return_invariant -- the refinement function provided
+    cells_to_refine_by -- an array, contents ignored
+
+    """
+    PS.depth += 1
+    PS_clear(PS)
+    cells_to_refine_by[0] = PS_split_point(PS, v)
+    return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
+
+cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, object S):
+    return 0
+
+cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
+    return 0
+
+cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S):
+    return 0
+
+def test_traverse_tree_trivially(int n=6, partition=[[0,1,2],[3,4],[5]],
+    canonical_label=True, base=False):
+    """
+    sage: tttt = sage.groups.perm_gps.partn_ref.tree_traversal.test_traverse_tree_trivially
+    sage: tttt()
+    12
+    sage: tttt(canonical_label=False, base=False)
+    12
+    sage: tttt(canonical_label=False, base=True)
+    12
+    sage: tttt(canonical_label=True, base=True)
+    12
+    sage: tttt(n=0, partition=[])
+    1
+    sage: tttt(n=0, partition=[], canonical_label=False, base=False)
+    1
+    sage: tttt(n=0, partition=[], canonical_label=False, base=True)
+    1
+    sage: tttt(n=0, partition=[], canonical_label=True, base=True)
+    1
+    
+    """
+    cdef int i, j
+    cdef traverse_return *output
+    cdef Integer I = Integer(0)
+    cdef int **part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
+    for i from 0 <= i < len(partition):
+        part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
+        for j from 0 <= j < len(partition[i]):
+            part[i][j] = partition[i][j]
+        part[i][len(partition[i])] = -1
+    part[len(partition)] = NULL
+    output = traverse_tree([], part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_structures_trivial, canonical_label, base, True)
+    mpz_set(I.value, output.order)
+    print I
+    for i from 0 <= i < len(partition):
+        sage_free(part[i])
+    sage_free(part)
+    mpz_clear(output.order)
+    sage_free(output.generators)
+    if base:
+        sage_free(output.base)
+    if canonical_label:
+        sage_free(output.relabeling)
+    sage_free(output)
+
+cdef traverse_return *traverse_tree(object S, int **partition, int n,
+    bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
+    int (*refine_and_return_invariant)\
+         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
+    int (*compare_structures)(int *gamma_1, int *gamma_2, object S),
+    bint canonical_label, bint base, bint order):
+    """
+    Traverse the search space for subgroup/canonical label calculation.
+    
+    INPUT:
+    S -- pointer to the structure
+    partition -- array representing a partition of the points
+    len_partition -- length of the partition
+    n -- the number of points (points are assumed to be 0,1,...,n-1)
+    canonical_label -- whether to search for canonical label - if True, return
+        the permutation taking S to its canonical label
+    base -- if True, return the base for which the given generating set is a
+        strong generating set
+    order -- if True, return the order of the automorphism group
+    all_children_are_equivalent -- pointer to a function
+        INPUT:
+        PS -- pointer to a partition stack
+        S -- pointer to the structure
+        OUTPUT:
+        bint -- returns True if it can be determined that all refinements below
+            the current one will result in an equivalent discrete partition
+    refine_and_return_invariant -- pointer to a function
+        INPUT:
+        PS -- pointer to a partition stack
+        S -- pointer to the structure
+        alpha -- an array consisting of numbers, which indicate the starting
+            positions of the cells to refine against (will likely be modified)
+        OUTPUT:
+        int -- returns an invariant under application of arbitrary permutations
+    compare_structures -- pointer to a function
+        INPUT:
+        gamma_1, gamma_2 -- (list) permutations of the points of S
+        S -- pointer to the structure
+        OUTPUT:
+        int -- 0 if gamma_1(S) = gamma_2(S), otherwise -1 or 1 (see docs for cmp),
+            such that the set of all structures is well-ordered
+    OUTPUT:
+    pointer to a traverse_return struct
+
+    """
+    cdef PartitionStack *current_ps, *first_ps, *label_ps
+    cdef int first_meets_current = -1
+    cdef int label_meets_current
+    cdef int current_kids_are_same = 1
+    cdef int first_kids_are_same
+    
+    cdef int *current_indicators, *first_indicators, *label_indicators
+    cdef int first_and_current_indicator_same
+    cdef int label_and_current_indicator_same = -1
+    cdef int compared_current_and_label_indicators
+    
+    cdef OrbitPartition *orbits_of_subgroup
+    cdef mpz_t subgroup_size
+    cdef int subgroup_primary_orbit_size = 0
+    cdef int minimal_in_primary_orbit
+    cdef int size_of_generator_array = 2*n*n
+    
+    cdef bitset_t *fixed_points_of_generators # i.e. fp
+    cdef bitset_t *minimal_cell_reps_of_generators # i.e. mcr
+    cdef int len_of_fp_and_mcr = 100
+    cdef int index_in_fp_and_mcr = -1
+    
+    cdef bitset_t *vertices_to_split
+    cdef bitset_t vertices_have_been_reduced
+    cdef int *permutation, *id_perm, *cells_to_refine_by
+    cdef int *vertices_determining_current_stack
+    
+    cdef int i, j, k
+    cdef bint discrete, automorphism, update_label
+    cdef bint backtrack, new_vertex, narrow, mem_err = 0
+    
+    cdef traverse_return *output = <traverse_return *> sage_malloc(sizeof(traverse_return))
+    if output is NULL:
+        raise MemoryError
+    
+    if n == 0:
+        output.generators = NULL
+        output.num_gens = 0
+        output.relabeling = NULL
+        output.base = NULL
+        output.base_size = 0
+        mpz_init_set_si(output.order, 1)
+        return output
+    
+    # Allocate:
+    mpz_init_set_si(subgroup_size, 1)
+    
+    output.generators = <int *> sage_malloc( size_of_generator_array * sizeof(int) )
+    output.num_gens = 0
+    if base:
+        output.base = <int *> sage_malloc( n * sizeof(int) )
+    
+    current_indicators = <int *> sage_malloc( n * sizeof(int) )
+    first_indicators = <int *> sage_malloc( n * sizeof(int) )
+    
+    if canonical_label:
+        label_indicators = <int *> sage_malloc( n * sizeof(int) )
+        output.relabeling = <int *> sage_malloc( n * sizeof(int) )
+
+    fixed_points_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
+    minimal_cell_reps_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
+    
+    vertices_to_split = <bitset_t *> sage_malloc( n * sizeof(bitset_t) )
+    permutation = <int *> sage_malloc( n * sizeof(int) )
+    id_perm = <int *> sage_malloc( n * sizeof(int) )
+    cells_to_refine_by = <int *> sage_malloc( n * sizeof(int) )
+    vertices_determining_current_stack = <int *> sage_malloc( n * sizeof(int) )
+    
+    current_ps = PS_new(n, False)
+    orbits_of_subgroup = OP_new(n)
+    
+    # Check for allocation failures:
+    cdef bint succeeded = 1
+    if canonical_label:
+        if label_indicators is NULL or output.relabeling is NULL:
+            succeeded = 0
+            if label_indicators is not NULL:
+                sage_free(label_indicators)
+            if output.relabeling is not NULL:
+                sage_free(output.relabeling)
+    if base:
+        if output.base is NULL:
+            succeeded = 0
+    if not succeeded                              or \
+       output.generators                  is NULL or \
+       current_indicators                 is NULL or \
+       first_indicators                   is NULL or \
+       fixed_points_of_generators         is NULL or \
+       minimal_cell_reps_of_generators    is NULL or \
+       vertices_to_split                  is NULL or \
+       permutation                        is NULL or \
+       id_perm                            is NULL or \
+       cells_to_refine_by                 is NULL or \
+       vertices_determining_current_stack is NULL or \
+       current_ps                         is NULL or \
+       orbits_of_subgroup                 is NULL:
+        if output.generators                  is not NULL:
+            sage_free(output.generators)
+        if base:
+            if output.base                        is not NULL:
+                sage_free(output.base)
+        if current_indicators                 is not NULL:
+            sage_free(current_indicators)
+        if first_indicators                   is not NULL:
+            sage_free(first_indicators)
+        if fixed_points_of_generators         is not NULL:
+            sage_free(fixed_points_of_generators)
+        if minimal_cell_reps_of_generators    is not NULL:
+            sage_free(minimal_cell_reps_of_generators)
+        if vertices_to_split                  is not NULL:
+            sage_free(vertices_to_split)
+        if permutation                        is not NULL:
+            sage_free(permutation)
+        if id_perm                            is not NULL:
+            sage_free(id_perm)
+        if cells_to_refine_by                 is not NULL:
+            sage_free(cells_to_refine_by)
+        if vertices_determining_current_stack is not NULL:
+            sage_free(vertices_determining_current_stack)
+        if current_ps is not NULL:
+            PS_dealloc(current_ps)
+        if orbits_of_subgroup is not NULL:
+            OP_dealloc(orbits_of_subgroup)
+        sage_free(output)
+        mpz_clear(subgroup_size)
+        raise MemoryError
+    
+    # Initialize bitsets, checking for allocation failures:
+    succeeded = 1
+    for i from 0 <= i < len_of_fp_and_mcr:
+        try:
+            bitset_init(fixed_points_of_generators[i], n)
+        except MemoryError:
+            succeeded = 0
+            for j from 0 <= j < i:
+                bitset_clear(fixed_points_of_generators[j])
+                bitset_clear(minimal_cell_reps_of_generators[j])
+            break
+        try:
+            bitset_init(minimal_cell_reps_of_generators[i], n)
+        except MemoryError:
+            succeeded = 0
+            for j from 0 <= j < i:
+                bitset_clear(fixed_points_of_generators[j])
+                bitset_clear(minimal_cell_reps_of_generators[j])
+            bitset_clear(fixed_points_of_generators[i])
+            break
+    if succeeded:
+        for i from 0 <= i < n:
+            try:
+                bitset_init(vertices_to_split[i], n)
+            except MemoryError:
+                succeeded = 0
+                for j from 0 <= j < i:
+                    bitset_clear(vertices_to_split[j])
+                for j from 0 <= j < len_of_fp_and_mcr:
+                    bitset_clear(fixed_points_of_generators[j])
+                    bitset_clear(minimal_cell_reps_of_generators[j])
+                break
+    if succeeded:
+        try:
+            bitset_init(vertices_have_been_reduced, n)
+        except MemoryError:
+            succeeded = 0
+            for j from 0 <= j < n:
+                bitset_clear(vertices_to_split[j])
+            for j from 0 <= j < len_of_fp_and_mcr:
+                bitset_clear(fixed_points_of_generators[j])
+                bitset_clear(minimal_cell_reps_of_generators[j])
+    if not succeeded:
+        sage_free(output.generators)
+        if base:
+            sage_free(output.base)
+        sage_free(current_indicators)
+        sage_free(first_indicators)
+        if canonical_label:
+            sage_free(output.relabeling)
+            sage_free(label_indicators)
+        sage_free(fixed_points_of_generators)
+        sage_free(minimal_cell_reps_of_generators)
+        sage_free(vertices_to_split)
+        sage_free(permutation)
+        sage_free(id_perm)
+        sage_free(cells_to_refine_by)
+        sage_free(vertices_determining_current_stack)
+        PS_dealloc(current_ps)
+        sage_free(output)
+        mpz_clear(subgroup_size)
+        raise MemoryError
+    
+    # Copy data of partition to current_ps.
+    i = 0
+    j = 0
+    while partition[i] is not NULL:
+        k = 0
+        while partition[i][k] != -1:
+            current_ps.entries[j+k] = partition[i][k]
+            current_ps.levels[j+k] = n
+            k += 1
+        current_ps.levels[j+k-1] = 0
+        PS_move_min_to_front(current_ps, j, j+k-1)
+        j += k
+        i += 1
+    current_ps.levels[j-1] = -1
+    current_ps.depth = 0
+    current_ps.degree = n
+    
+    # default values of "infinity"
+    for i from 0 <= i < n:
+        first_indicators[i] = -1
+    if canonical_label:
+        for i from 0 <= i < n:
+            label_indicators[i] = -1
+    
+    # set up the identity permutation
+    for i from 0 <= i < n:
+        id_perm[i] = i
+    
+    # Our first refinement needs to check every cell of the partition,
+    # so cells_to_refine_by needs to be a list of pointers to each cell.
+    j = 1
+    cells_to_refine_by[0] = 0
+    for i from 0 < i < n:
+        if current_ps.levels[i-1] == 0:
+            cells_to_refine_by[j] = i
+            j += 1
+    # Ignore the invariant this time, since we are
+    # creating the root of the search tree.
+    refine_and_return_invariant(current_ps, S, cells_to_refine_by, j)
+
+    # Refine down to a discrete partition
+    while not PS_is_discrete(current_ps):
+        if not all_children_are_equivalent(current_ps, S):
+            current_kids_are_same = current_ps.depth + 1
+        vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
+        bitset_unset(vertices_have_been_reduced, current_ps.depth)
+        i = current_ps.depth
+        current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
+        first_indicators[i] = current_indicators[i]
+        if canonical_label:
+            label_indicators[i] = current_indicators[i]
+        if base:
+            output.base[i] = vertices_determining_current_stack[i]
+
+    first_meets_current = current_ps.depth
+    first_kids_are_same = current_ps.depth
+    first_and_current_indicator_same = current_ps.depth
+    first_ps = PS_copy(current_ps)
+    if canonical_label:
+        label_meets_current = current_ps.depth
+        label_and_current_indicator_same = current_ps.depth
+        compared_current_and_label_indicators = 0
+        label_ps = PS_copy(current_ps)
+    current_ps.depth -= 1
+    
+    # Main loop:
+    while current_ps.depth != -1:
+        
+        # If necessary, update label_meets_current
+        if canonical_label and label_meets_current > current_ps.depth:
+            label_meets_current = current_ps.depth
+
+        # I. Search for a new vertex to split, and update subgroup information
+
+        new_vertex = 0
+        if current_ps.depth > first_meets_current:
+            # If we are not at a node of the first stack, reduce size of
+            # vertices_to_split by using the symmetries we already know.
+            if not bitset_check(vertices_have_been_reduced, current_ps.depth):
+                for i from 0 <= i <= index_in_fp_and_mcr:
+                    j = 0
+                    while j < current_ps.depth and bitset_check(fixed_points_of_generators[i], vertices_determining_current_stack[j]):
+                        j += 1
+                    # If each vertex split so far is fixed by generator i,
+                    # then remove elements of vertices_to_split which are
+                    # not minimal in their orbits under generator i.
+                    if j == current_ps.depth:
+                        for k from 0 <= k < n:
+                            if bitset_check(vertices_to_split[current_ps.depth], k) and not bitset_check(minimal_cell_reps_of_generators[i], k):
+                                bitset_flip(vertices_to_split[current_ps.depth], k)
+                bitset_flip(vertices_have_been_reduced, current_ps.depth)
+            # Look for a new point to split.
+            i = vertices_determining_current_stack[current_ps.depth] + 1
+            while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
+                i += 1
+            if i < n:
+                # There is a new point.
+                vertices_determining_current_stack[current_ps.depth] = i
+                new_vertex = 1
+            else:
+                # No new point: backtrack.
+                current_ps.depth -= 1
+        else:
+            # If we are at a node of the first stack, the above reduction
+            # will not help. Also, we must update information about
+            # primary orbits here.
+            if current_ps.depth < first_meets_current:
+                # If we are done searching under this part of the first
+                # stack, then first_meets_current is one higher, and we
+                # are looking at a new primary orbit (corresponding to a
+                # larger subgroup in the stabilizer chain).
+                first_meets_current = current_ps.depth
+                for i from 0 <= i < n:
+                    if bitset_check(vertices_to_split[current_ps.depth], i):
+                        minimal_in_primary_orbit = i
+                        break
+            while 1:
+                i = vertices_determining_current_stack[current_ps.depth]
+                # This was the last point to be split here.
+                # If it is in the same orbit as minimal_in_primary_orbit,
+                # then count it as an element of the primary orbit.
+                if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit):
+                    subgroup_primary_orbit_size += 1
+                # Look for a new point to split.
+                i += 1
+                while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
+                    i += 1
+                if i < n:
+                    # There is a new point.
+                    vertices_determining_current_stack[current_ps.depth] = i
+                    if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i:
+                        new_vertex = 1
+                        break
+                else:
+                    # No new point: backtrack.
+                    # Note that now, we are backtracking up the first stack.
+                    vertices_determining_current_stack[current_ps.depth] = -1
+                    # If every choice of point to split gave something in the
+                    # (same!) primary orbit, then all children of the first
+                    # stack at this point are equivalent.
+                    j = 0
+                    for i from 0 <= i < n:
+                        if bitset_check(vertices_to_split[current_ps.depth], i):
+                            j += 1
+                    if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1:
+                        first_kids_are_same = current_ps.depth
+                    # Update group size and backtrack.
+                    mpz_mul_si(subgroup_size, subgroup_size, subgroup_primary_orbit_size)
+                    subgroup_primary_orbit_size = 0
+                    current_ps.depth -= 1
+                    break
+        if not new_vertex:
+            continue
+
+        if current_kids_are_same > current_ps.depth + 1:
+            current_kids_are_same = current_ps.depth + 1
+        if first_and_current_indicator_same > current_ps.depth:
+            first_and_current_indicator_same = current_ps.depth
+        if canonical_label and label_and_current_indicator_same > current_ps.depth:
+            label_and_current_indicator_same = current_ps.depth
+            compared_current_and_label_indicators = 0
+        
+        # II. Refine down to a discrete partition, or until
+        # we leave the part of the tree we are interested in
+        discrete = 0
+        while 1:
+            i = current_ps.depth
+            current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
+            if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]:
+                first_and_current_indicator_same += 1
+            if canonical_label:
+                if compared_current_and_label_indicators == 0:
+                    if label_indicators[i] == -1:
+                        compared_current_and_label_indicators = -1
+                    else:
+                        compared_current_and_label_indicators = cmp(current_indicators[i], label_indicators[i])
+                if label_and_current_indicator_same == i and compared_current_and_label_indicators == 0:
+                    label_and_current_indicator_same += 1
+                if compared_current_and_label_indicators > 0:
+                    label_indicators[i] = current_indicators[i]
+            if first_and_current_indicator_same < current_ps.depth and (not canonical_label or compared_current_and_label_indicators < 0):
+                break
+            if PS_is_discrete(current_ps):
+                discrete = 1
+                break
+            vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
+            bitset_unset(vertices_have_been_reduced, current_ps.depth)
+            if not all_children_are_equivalent(current_ps, S):
+                current_kids_are_same = current_ps.depth + 1
+
+        # III. Check for automorphisms and labels
+        automorphism = 0
+        backtrack = 1 - discrete
+        if discrete:
+            if current_ps.depth == first_and_current_indicator_same:
+                PS_get_perm_from(current_ps, first_ps, permutation)
+                if compare_structures(permutation, id_perm, S) == 0:
+                    automorphism = 1
+        if not automorphism:
+            if (not canonical_label) or (compared_current_and_label_indicators < 0):
+                backtrack = 1
+            else:
+                # if we get here, discrete must be true
+                update_label = 0
+                if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth):
+                    update_label = 1
+                else:
+                    i = compare_structures(current_ps.entries, label_ps.entries, S)
+                    if i > 0:
+                        update_label = 1
+                    elif i < 0:
+                        backtrack = 1
+                    else:
+                        automorphism = 1
+                if update_label:
+                    PS_copy_from_to(current_ps, label_ps)
+                    compared_current_and_label_indicators = 0
+                    label_meets_current = current_ps.depth
+                    label_and_current_indicator_same = current_ps.depth
+                    label_indicators[current_ps.depth] = -1
+                    backtrack = 1
+        if automorphism:
+            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
+                index_in_fp_and_mcr += 1
+            for i from 0 <= i < n:
+                if permutation[i] == i:
+                    bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
+                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+                else:
+                    bitset_unset(fixed_points_of_generators[index_in_fp_and_mcr], i)
+                    k = i
+                    j = permutation[i]
+                    while j != i:
+                        if j < k: k = j
+                        j = permutation[j]
+                    if k == i:
+                        bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+                    else:
+                        bitset_unset(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+            if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser
+                if n*output.num_gens == size_of_generator_array:
+                    # must double its size
+                    size_of_generator_array *= 2
+                    output.generators = <int *> sage_realloc( output.generators, size_of_generator_array )
+                    if output.generators is NULL:
+                        mem_err = True
+                        break # main loop
+                j = n*output.num_gens
+                for i from 0 <= i < n:
+                    output.generators[j + i] = permutation[i]
+                output.num_gens += 1
+                if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit:
+                    current_ps.depth = first_meets_current
+                    continue # main loop
+            if canonical_label:
+                current_ps.depth = label_meets_current
+            else:
+                current_ps.depth = first_meets_current
+            if bitset_check(vertices_have_been_reduced, current_ps.depth):
+                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+            continue # main loop
+        if backtrack:
+            i = current_ps.depth
+            current_ps.depth = min(max(first_kids_are_same - 1, label_and_current_indicator_same), current_kids_are_same - 1)
+            if i == current_kids_are_same:
+                continue # main loop
+            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
+                index_in_fp_and_mcr += 1
+            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
+            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+            j = current_ps.depth
+            current_ps.depth = i # just for mcr and fixed functions...
+            for i from 0 <= i < n:
+                if PS_is_mcr(current_ps, i):
+                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
+                    if PS_is_fixed(current_ps, i):
+                        bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
+            current_ps.depth = j
+            if bitset_check(vertices_have_been_reduced, current_ps.depth):
+                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
+
+    # End of main loop.
+
+    mpz_init_set(output.order, subgroup_size)
+    if canonical_label:
+        for i from 0 <= i < n:
+            output.relabeling[label_ps.entries[i]] = i
+    
+    # Deallocate:
+    for i from 0 <= i < len_of_fp_and_mcr:
+        bitset_clear(fixed_points_of_generators[i])
+        bitset_clear(minimal_cell_reps_of_generators[i])
+    for i from 0 <= i < n:
+        bitset_clear(vertices_to_split[i])
+    bitset_clear(vertices_have_been_reduced)
+    sage_free(current_indicators)
+    sage_free(first_indicators)
+    if canonical_label:
+        PS_dealloc(label_ps)
+        sage_free(label_indicators)
+    sage_free(fixed_points_of_generators)
+    sage_free(minimal_cell_reps_of_generators)
+    sage_free(vertices_to_split)
+    sage_free(permutation)
+    sage_free(id_perm)
+    sage_free(cells_to_refine_by)
+    sage_free(vertices_determining_current_stack)
+    PS_dealloc(current_ps)
+    PS_dealloc(first_ps)
+    OP_dealloc(orbits_of_subgroup)
+    mpz_clear(subgroup_size)
+    if not mem_err:
+        return output
+    else:
+        mpz_clear(output.order)
+        sage_free(output.generators)
+        if base:
+            sage_free(output.base)
+        if canonical_label:
+            sage_free(output.relabeling)
+        sage_free(output)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

setup.py

@@ -675 +675 @@
 
     Extension('sage.groups.perm_gps.permgroup_element',
               sources = ['sage/groups/perm_gps/permgroup_element.pyx']), \
+
+    Extension('sage.groups.perm_gps.partn_ref.tree_traversal',
+              sources = ['sage/groups/perm_gps/partn_ref/tree_traversal.pyx']), \
+
+    Extension('sage.groups.perm_gps.partn_ref.refinement_graphs',
+              sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx']), \
 
     Extension('sage.structure.sage_object',
               sources = ['sage/structure/sage_object.pyx']), \
@@ -1370 +1376 @@
                      'sage.groups.abelian_gps',
                      'sage.groups.matrix_gps',
                      'sage.groups.perm_gps',
+                     'sage.groups.perm_gps.partn_ref',
                      
                      'sage.interfaces',