Ticket #3676: trac3676-refactor_graph_isom.patch

File trac3676-refactor_graph_isom.patch, 79.4 KB (added by rlm, 13 years ago)
  • sage/graphs/graph_generators.py

    # HG changeset patch
    # User Robert L. Miller <rlm@rlmiller.org>
    # Date 1216425781 25200
    # Node ID 074731a2698613667ad1b649a96a8446a892ffd6
    # Parent  91af4c3f6b92316e4117af342d82be33cf5d949b
    Refactored graph isomorphism code
    
    diff -r 91af4c3f6b92 -r 074731a26986 sage/graphs/graph_generators.py
    a b class DiGraphGenerators(): 
    29872987        import networkx
    29882988        return graph.DiGraph(networkx.gnc_graph(n, seed))
    29892989
     2990    def RandomDirectedGNP(self, n, p):
     2991        r"""
     2992        Returns a random digraph on $n$ nodes.  Each edge is inserted
     2993        independently with probability $p$.
     2994
     2995        REFERENCES:
     2996            [1] P. Erdos and A. Renyi, On Random Graphs, Publ. Math. 6, 290 (1959).
     2997            [2] E. N. Gilbert, Random Graphs, Ann. Math. Stat., 30, 1141 (1959).
     2998       
     2999        PLOTTING:
     3000        When plotting, this graph will use the default spring-layout
     3001        algorithm, unless a position dictionary is specified.
     3002               
     3003        EXAMPLE:
     3004            sage: digraphs.RandomDirectedGNP(10, .2).num_verts()
     3005            10
     3006       
     3007        """
     3008        from random import random
     3009        D = graph.DiGraph(n)
     3010        for i in xrange(n):
     3011            for j in xrange(i):
     3012                if random() < p:
     3013                    D.add_edge(i,j)
     3014            for j in xrange(i+1,n):
     3015                if random() < p:
     3016                    D.add_edge(i,j)
     3017        return D
     3018
    29903019    def RandomDirectedGNR(self, n, p, seed=None):
    29913020        """
    29923021        Returns a random GNR (growing network with redirection) digraph with n
  • new file sage/groups/perm_gps/partn_ref/__init__.py

    diff -r 91af4c3f6b92 -r 074731a26986 sage/groups/perm_gps/partn_ref/__init__.py
    - +  
     1
  • new file sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

    diff -r 91af4c3f6b92 -r 074731a26986 sage/groups/perm_gps/partn_ref/refinement_graphs.pxd
    - +  
     1
     2#*****************************************************************************
     3#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     4#
     5# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     6#                         http://www.gnu.org/licenses/
     7#*****************************************************************************
     8
     9include '../../../ext/cdefs.pxi'
     10include '../../../ext/stdsage.pxi'
     11
     12from sage.graphs.base.c_graph cimport CGraph
     13from sage.graphs.base.sparse_graph cimport SparseGraph
     14from sage.graphs.base.dense_graph cimport DenseGraph
     15from sage.rings.integer cimport Integer
     16from 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
     17
     18cdef class GraphStruct:
     19    cdef CGraph G
     20    cdef bint directed
     21    cdef bint use_indicator
     22    cdef int *scratch # length 3n+1
     23
     24cdef int refine_by_degree(PartitionStack *, object, int *, int)
     25cdef int compare_graphs(int *, int *, object)
     26cdef bint all_children_are_equivalent(PartitionStack *, object)
     27cdef int degree(PartitionStack *, CGraph, int, int, bint)
     28cdef int sort_by_function(PartitionStack *, int, int *)
     29
     30
     31
  • new file sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

    diff -r 91af4c3f6b92 -r 074731a26986 sage/groups/perm_gps/partn_ref/refinement_graphs.pyx
    - +  
     1"""
     2Module for graph-theoretic partition backtrack functions.
     3
     4DOCTEST:
     5    sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
     6
     7"""
     8
     9#*****************************************************************************
     10#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     11#
     12# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     13#                         http://www.gnu.org/licenses/
     14#*****************************************************************************
     15
     16def search_tree(G_in, partition, lab=True, dig=False, dict_rep=False, certify=False,
     17                    verbosity=0, use_indicator_function=True, sparse=False,
     18                    base=False, order=False):
     19    """
     20    Compute canonical labels and automorphism groups of graphs.
     21   
     22    INPUT:
     23    G_in -- a Sage graph
     24    partition -- a list of lists representing a partition of the vertices
     25    lab -- if True, return the canonical label in addition to the automorphism
     26        group.
     27    dig -- if True, does not use Lemma 2.25 in [1], and the algorithm is valid
     28        for digraphs and graphs with loops.
     29    dict_rep -- if True, explain which vertices are which elements of the set
     30        {1,2,...,n} in the representation of the automorphism group.
     31    certify -- if True, return the relabeling from G to its canonical
     32        label.
     33    verbosity -- currently ignored
     34    use_indicator_function -- option to turn off indicator function
     35        (False -> slower)
     36    sparse -- whether to use sparse or dense representation of the graph
     37        (ignored if G is already a CGraph - see sage.graphs.base)
     38    base -- whether to return the first sequence of split vertices (used in
     39        computing the order of the group)
     40    order -- whether to return the order of the automorphism group
     41
     42    OUTPUT:
     43    Depends on the options. If more than one thing is returned, they are in a
     44    tuple in the following order:
     45        list of generators in list-permutation format -- always
     46        dict -- if dict_rep
     47        graph -- if lab
     48        dict -- if certify
     49        list -- if base
     50        integer -- if order
     51
     52    DOCTEST:
     53        sage: st = sage.groups.perm_gps.partn_ref.refinement_graphs.search_tree
     54        sage: from sage.graphs.base.dense_graph import DenseGraph
     55        sage: from sage.graphs.base.sparse_graph import SparseGraph
     56
     57    Graphs on zero vertices:
     58        sage: G = Graph()
     59        sage: st(G, [[]], order=True)
     60        ([], Graph on 0 vertices, 1)
     61
     62    Graphs on one vertex:
     63        sage: G = Graph(1)
     64        sage: st(G, [[0]], order=True)
     65        ([], Graph on 1 vertex, 1)
     66
     67    Graphs on two vertices:
     68        sage: G = Graph(2)
     69        sage: st(G, [[0,1]], order=True)
     70        ([[1, 0]], Graph on 2 vertices, 2)
     71        sage: st(G, [[0],[1]], order=True)
     72        ([], Graph on 2 vertices, 1)
     73        sage: G.add_edge(0,1)
     74        sage: st(G, [[0,1]], order=True)
     75        ([[1, 0]], Graph on 2 vertices, 2)
     76        sage: st(G, [[0],[1]], order=True)
     77        ([], Graph on 2 vertices, 1)
     78
     79    Graphs on three vertices:
     80        sage: G = Graph(3)
     81        sage: st(G, [[0,1,2]], order=True)
     82        ([[0, 2, 1], [1, 0, 2]], Graph on 3 vertices, 6)
     83        sage: st(G, [[0],[1,2]], order=True)
     84        ([[0, 2, 1]], Graph on 3 vertices, 2)
     85        sage: st(G, [[0],[1],[2]], order=True)
     86        ([], Graph on 3 vertices, 1)
     87        sage: G.add_edge(0,1)
     88        sage: st(G, [range(3)], order=True)
     89        ([[1, 0, 2]], Graph on 3 vertices, 2)
     90        sage: st(G, [[0],[1,2]], order=True)
     91        ([], Graph on 3 vertices, 1)
     92        sage: st(G, [[0,1],[2]], order=True)
     93        ([[1, 0, 2]], Graph on 3 vertices, 2)
     94
     95    The Dodecahedron has automorphism group of size 120:
     96        sage: G = graphs.DodecahedralGraph()
     97        sage: Pi = [range(20)]
     98        sage: st(G, Pi, order=True)[2]
     99        120
     100
     101    The three-cube has automorphism group of size 48:
     102        sage: G = graphs.CubeGraph(3)
     103        sage: G.relabel()
     104        sage: Pi = [G.vertices()]
     105        sage: st(G, Pi, order=True)[2]
     106        48
     107
     108    We obtain the same output using different types of Sage graphs:
     109        sage: G = graphs.DodecahedralGraph()
     110        sage: GD = DenseGraph(20)
     111        sage: GS = SparseGraph(20)
     112        sage: for i,j,_ in G.edge_iterator():
     113        ...    GD.add_arc(i,j); GD.add_arc(j,i)
     114        ...    GS.add_arc(i,j); GS.add_arc(j,i)
     115        sage: Pi=[range(20)]
     116        sage: a,b = st(G, Pi)
     117        sage: asp,bsp = st(GS, Pi)
     118        sage: ade,bde = st(GD, Pi)
     119        sage: bsg = Graph(implementation='networkx')
     120        sage: bdg = Graph(implementation='networkx')
     121        sage: for i in range(20):
     122        ...    for j in range(20):
     123        ...        if bsp.has_arc(i,j):
     124        ...            bsg.add_edge(i,j)
     125        ...        if bde.has_arc(i,j):
     126        ...            bdg.add_edge(i,j)
     127        sage: print a, b.graph6_string()
     128        [[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
     129        sage: a == asp
     130        True
     131        sage: a == ade
     132        True
     133        sage: b == bsg
     134        True
     135        sage: b == bdg
     136        True
     137
     138    Cubes!
     139        sage: C = graphs.CubeGraph(1)
     140        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
     141        2
     142        sage: C = graphs.CubeGraph(2)
     143        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
     144        8
     145        sage: C = graphs.CubeGraph(3)
     146        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
     147        48
     148        sage: C = graphs.CubeGraph(4)
     149        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
     150        384
     151        sage: C = graphs.CubeGraph(5)
     152        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
     153        3840
     154        sage: C = graphs.CubeGraph(6)
     155        sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order
     156        46080
     157
     158    One can also turn off the indicator function (note- this will take longer)
     159        sage: D1 = DiGraph({0:[2],2:[0],1:[1]}, loops=True)
     160        sage: D2 = DiGraph({1:[2],2:[1],0:[0]}, loops=True)
     161        sage: a,b = st(D1, [D1.vertices()], dig=True, use_indicator_function=False)
     162        sage: c,d = st(D2, [D2.vertices()], dig=True, use_indicator_function=False)
     163        sage: b==d
     164        True
     165
     166    This example is due to Chris Godsil:
     167        sage: HS = graphs.HoffmanSingletonGraph()
     168        sage: clqs = (HS.complement()).cliques()
     169        sage: alqs = [Set(c) for c in clqs if len(c) == 15]
     170        sage: Y = Graph([alqs, lambda s,t: len(s.intersection(t))==0], implementation='networkx')
     171        sage: Y0,Y1 = Y.connected_components_subgraphs()
     172        sage: st(Y0, [Y0.vertices()])[1] == st(Y1, [Y1.vertices()])[1]
     173        True
     174        sage: st(Y0, [Y0.vertices()])[1] == st(HS, [HS.vertices()])[1]
     175        True
     176        sage: st(HS, [HS.vertices()])[1] == st(Y1, [Y1.vertices()])[1]
     177        True
     178
     179    Certain border cases need to be tested as well:
     180        sage: G = Graph('Fll^G')
     181        sage: a,b,c = st(G, [range(G.num_verts())], order=True); b
     182        Graph on 7 vertices
     183        sage: c
     184        48
     185        sage: G = Graph(21)
     186        sage: st(G, [range(G.num_verts())], order=True)[2] == factorial(21)
     187        True
     188
     189        sage: G = Graph('^????????????????????{??N??@w??FaGa?PCO@CP?AGa?_QO?Q@G?CcA??cc????Bo????{????F_')
     190        sage: perm = {3:15, 15:3}
     191        sage: H = G.relabel(perm, inplace=False)
     192        sage: st(G, [range(G.num_verts())])[1] == st(H, [range(H.num_verts())])[1]
     193        True
     194
     195    """
     196    cdef CGraph G
     197    cdef int i, j, n
     198    cdef Integer I
     199    cdef traverse_return *output
     200    cdef int **part
     201    from sage.graphs.graph import GenericGraph, Graph, DiGraph
     202    if isinstance(G_in, GenericGraph):
     203        n = G_in.num_verts()
     204        G_in = G_in.copy()
     205        if G_in.vertices() != range(n):
     206            to = G_in.relabel(return_map=True)
     207            frm = {}
     208            for v in to.iterkeys():
     209                frm[to[v]] = v
     210            partition = [[to[v] for v in cell] for cell in partition]
     211        else:
     212            to = range(n)
     213            frm = to
     214        if sparse:
     215            G = SparseGraph(n)
     216        else:
     217            G = DenseGraph(n)
     218        if G_in.is_directed():
     219            for i from 0 <= i < n:
     220                for _,j,_ in G_in.outgoing_edge_iterator(i):
     221                    G.add_arc(i,j)
     222        else:
     223            for i from 0 <= i < n:
     224                for _,j,_ in G_in.edge_iterator(i):
     225                    if j <= i:
     226                        G.add_arc(i,j)
     227                        G.add_arc(j,i)
     228    elif isinstance(G_in, CGraph):
     229        G = <CGraph> G_in
     230        n = G.num_verts
     231        to = range(n)
     232        frm = to
     233    else:
     234        raise TypeError("G must be a Sage graph.")
     235   
     236    part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     237    if part is NULL:
     238        raise MemoryError
     239    for i from 0 <= i < len(partition):
     240        part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
     241        if part[i] is NULL:
     242            raise MemoryError
     243        for j from 0 <= j < len(partition[i]):
     244            part[i][j] = partition[i][j]
     245        part[i][len(partition[i])] = -1
     246    part[len(partition)] = NULL
     247
     248    cdef GraphStruct GS = GraphStruct()
     249    GS.G = G
     250    GS.directed = 1 if dig else 0
     251    GS.use_indicator = 1 if use_indicator_function else 0
     252    GS.scratch = <int *> sage_malloc( (3*G.num_verts + 1) * sizeof(int) )
     253    if GS.scratch is NULL:
     254        raise MemoryError
     255
     256    output = traverse_tree(GS, part, G.num_verts, &all_children_are_equivalent, &refine_by_degree, &compare_graphs, lab, base, order)
     257    sage_free( GS.scratch )
     258    # prepare output
     259    list_of_gens = []
     260    for i from 0 <= i < output.num_gens:
     261        list_of_gens.append([output.generators[j+i*G.num_verts] for j from 0 <= j < G.num_verts])
     262    return_tuple = [list_of_gens]
     263    if dict_rep:
     264        ddd = {}
     265        for v in frm.iterkeys():
     266            if frm[v] != 0:
     267                ddd[v] = frm[v]
     268            else:
     269                ddd[v] = n
     270        return_tuple.append(ddd)
     271    if lab:
     272        if isinstance(G_in, GenericGraph):
     273            G_C = G_in.copy()
     274            G_C.relabel([output.relabeling[i] for i from 0 <= i < n])
     275        else:
     276            if isinstance(G, SparseGraph):
     277                G_C = SparseGraph(n)
     278            else:
     279                G_C = DenseGraph(n)
     280            for i from 0 <= i < n:
     281                for j in G.out_neighbors(i):
     282                    G_C.add_arc(output.relabeling[i],output.relabeling[j])
     283        return_tuple.append(G_C)
     284    if certify:
     285        dd = {}
     286        for i from 0 <= i < G.num_verts:
     287            dd[frm[i]] = output.relabeling[i]
     288        return_tuple.append(dd)
     289    if base:
     290        return_tuple.append([output.base[i] for i from 0 <= i < output.base_size])
     291    if order:
     292        I = Integer()
     293        mpz_set(I.value, output.order)
     294        return_tuple.append(I)
     295    for i from 0 <= i < len(partition):
     296        sage_free(part[i])
     297    sage_free(part)
     298    mpz_clear(output.order)
     299    sage_free(output.generators)
     300    if base:
     301        sage_free(output.base)
     302    if lab:
     303        sage_free(output.relabeling)
     304    sage_free(output)
     305    if len(return_tuple) == 1:
     306        return return_tuple[0]
     307    else:
     308        return tuple(return_tuple)
     309
     310cdef int refine_by_degree(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     311    """
     312    Refines the input partition by checking degrees of vertices to the given
     313    cells.
     314   
     315    INPUT:
     316    PS -- a partition stack, whose finest partition is the partition to be
     317        refined.
     318    S -- a graph struct object, which contains scratch space, the graph in
     319        question, and some flags.
     320    cells_to_refine_by -- a list of pointers to cells to check degrees against
     321        in refining the other cells (updated in place)
     322    ctrb_len -- how many cells in the above
     323   
     324    OUTPUT:
     325    An invariant, such that if the function is called with a permutation of the
     326    situation (i.e. there is a gamma taking one graph to the other, which also
     327    takes one partition to the other), the same number is returned.
     328    """
     329    cdef GraphStruct GS = <GraphStruct> S
     330    cdef CGraph G = GS.G
     331    cdef int current_cell_against = 0
     332    cdef int current_cell, i, r
     333    cdef int first_largest_subcell
     334    cdef int invariant = 1
     335    cdef int *degrees = GS.scratch # length 3n+1
     336    cdef bint necessary_to_split_cell
     337    cdef int against_index
     338    while not PS_is_discrete(PS) and current_cell_against < ctrb_len:
     339        invariant += 1
     340        current_cell = 0
     341        while current_cell < PS.degree:
     342            invariant += 50
     343            i = current_cell
     344            necessary_to_split_cell = 0
     345            while 1:
     346                degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 0)
     347                if degrees[i-current_cell] != degrees[0]:
     348                    necessary_to_split_cell = 1
     349                i += 1
     350                if PS.levels[i-1] <= PS.depth:
     351                    break
     352            # now, i points to the next cell (before refinement)
     353            if necessary_to_split_cell:
     354                invariant += 10
     355                first_largest_subcell = sort_by_function(PS, current_cell, degrees)
     356                invariant += first_largest_subcell
     357                against_index = current_cell_against
     358                while against_index < ctrb_len:
     359                    if cells_to_refine_by[against_index] == current_cell:
     360                        cells_to_refine_by[against_index] = first_largest_subcell
     361                        break
     362                    against_index += 1
     363                against_index = ctrb_len
     364                r = current_cell
     365                while 1:
     366                    if r == current_cell or PS.levels[r-1] == PS.depth:
     367                        if r != first_largest_subcell:
     368                            cells_to_refine_by[against_index] = r
     369                            against_index += 1
     370                            ctrb_len += 1
     371                    r += 1
     372                    if r >= i:
     373                        break
     374                invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 0)
     375                invariant += (i - current_cell)
     376            current_cell = i
     377        if GS.directed:
     378            # if we are looking at a digraph, also compute
     379            # the reverse degrees and sort by them
     380            current_cell = 0
     381            while current_cell < PS.degree: # current_cell is still a valid cell
     382                invariant += 20
     383                i = current_cell
     384                necessary_to_split_cell = 0
     385                while 1:
     386                    degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 1)
     387                    if degrees[i-current_cell] != degrees[0]:
     388                        necessary_to_split_cell = 1
     389                    i += 1
     390                    if PS.levels[i-1] <= PS.depth:
     391                        break
     392                # now, i points to the next cell (before refinement)
     393                if necessary_to_split_cell:
     394                    invariant += 7
     395                    first_largest_subcell = sort_by_function(PS, current_cell, degrees)
     396                    invariant += first_largest_subcell
     397                    against_index = current_cell_against
     398                    while against_index < ctrb_len:
     399                        if cells_to_refine_by[against_index] == current_cell:
     400                            cells_to_refine_by[against_index] = first_largest_subcell
     401                            break
     402                        against_index += 1
     403                    against_index = ctrb_len
     404                    r = current_cell
     405                    while 1:
     406                        if r == current_cell or PS.levels[r-1] == PS.depth:
     407                            if r != first_largest_subcell:
     408                                cells_to_refine_by[against_index] = r
     409                                against_index += 1
     410                                ctrb_len += 1
     411                        r += 1
     412                        if r >= i:
     413                            break
     414                    invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 1)
     415                    invariant += (i - current_cell)
     416                current_cell = i
     417        current_cell_against += 1
     418    if GS.use_indicator:
     419        return invariant
     420    else:
     421        return 0
     422   
     423cdef int compare_graphs(int *gamma_1, int *gamma_2, object S):
     424    """
     425    Compare gamma_1(G) and gamma_2(G), returning 0 iff gamma_1(G) == gamma_2(G).
     426   
     427    INPUT:
     428    gamma_1, gamma_2 -- list permutations (inverse)
     429    S -- a graph struct object
     430   
     431    """
     432    cdef int i, j
     433    cdef GraphStruct GS = <GraphStruct> S
     434    cdef CGraph G = GS.G
     435    for i from 0 <= i < G.num_verts:
     436        for j from 0 <= j < G.num_verts:
     437            if G.has_arc_unsafe(gamma_1[i], gamma_1[j]):
     438                if not G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
     439                    return 1
     440            elif G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
     441                return -1
     442    return 0
     443
     444cdef bint all_children_are_equivalent(PartitionStack *PS, object S):
     445    """
     446    Returns True if any refinement of the current partition results in the same
     447    structure. Converse does not hold in general.
     448   
     449    INPUT:
     450    PS -- the partition stack to be checked
     451    S -- a graph struct object
     452    """
     453    cdef GraphStruct GS = <GraphStruct> S
     454    cdef CGraph G = GS.G
     455    cdef int i, n = PS.degree
     456    cdef bint in_cell = 0
     457    cdef int nontrivial_cells = 0
     458    cdef int total_cells = PS_num_cells(PS)
     459    if n <= total_cells + 4:
     460        return 1
     461    for i from 0 <= i < n-1:
     462        if PS.levels[i] <= PS.depth:
     463            if in_cell:
     464                nontrivial_cells += 1
     465            in_cell = 0
     466        else:
     467            in_cell = 1
     468    if in_cell:
     469        nontrivial_cells += 1
     470    if n == total_cells + nontrivial_cells:
     471        return 1
     472    if n == total_cells + nontrivial_cells + 1:
     473        return 1
     474    return 0
     475
     476cdef int degree(PartitionStack *PS, CGraph G, int entry, int cell_index, bint reverse):
     477    """
     478    Returns the number of edges from the vertex corresponding to entry to
     479    vertices in the cell corresponding to cell_index.
     480   
     481    INPUT:
     482    PS -- the partition stack to be checked
     483    S -- a graph struct object
     484    entry -- the position of the vertex in question in the entries of PS
     485    cell_index -- the starting position of the cell in question in the entries
     486        of PS
     487    reverse -- whether to check for arcs in the other direction
     488    """
     489    cdef int num_arcs = 0
     490    entry = PS.entries[entry]
     491    if not reverse:
     492        while 1:
     493            if G.has_arc_unsafe(PS.entries[cell_index], entry):
     494                num_arcs += 1
     495            if PS.levels[cell_index] > PS.depth:
     496                cell_index += 1
     497            else:
     498                break
     499    else:
     500        while 1:
     501            if G.has_arc_unsafe(entry, PS.entries[cell_index]):
     502                num_arcs += 1
     503            if PS.levels[cell_index] > PS.depth:
     504                cell_index += 1
     505            else:
     506                break
     507    return num_arcs
     508
     509cdef int sort_by_function(PartitionStack *PS, int start, int *degrees):
     510    """
     511    A simple counting sort, given the degrees of vertices to a certain cell.
     512   
     513    INPUT:
     514    PS -- the partition stack to be checked
     515    start -- beginning index of the cell to be sorted
     516    degrees -- the values to be sorted by
     517
     518    """
     519    cdef int n = PS.degree
     520    cdef int i, j, m = 2*n, max, max_location
     521    cdef int *counts = degrees + n, *output = degrees + 2*n + 1
     522    for i from 0 <= i <= n:
     523        counts[i] = 0
     524    i = 0
     525    while PS.levels[i+start] > PS.depth:
     526        counts[degrees[i]] += 1
     527        i += 1
     528    counts[degrees[i]] += 1
     529    # i+start is the right endpoint of the cell now
     530    max = counts[0]
     531    max_location = 0
     532    for j from 0 < j <= n:
     533        if counts[j] > max:
     534            max = counts[j]
     535            max_location = j
     536        counts[j] += counts[j - 1]
     537    for j from i >= j >= 0:
     538        counts[degrees[j]] -= 1
     539        output[counts[degrees[j]]] = PS.entries[start+j]
     540    max_location = counts[max_location]+start
     541    for j from 0 <= j <= i:
     542        PS.entries[start+j] = output[j]
     543    j = 1
     544    while j <= n and counts[j] <= i:
     545        if counts[j] > 0:
     546            PS.levels[start + counts[j] - 1] = PS.depth
     547        PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1)
     548        j += 1
     549    return max_location
     550   
     551   
     552def all_labeled_graphs(n):
     553    """
     554    Returns all labeled graphs on n vertices {0,1,...,n-1}. Used in
     555    classifying isomorphism types (naive approach), and more importantly
     556    in benchmarking the search algorithm.
     557   
     558    EXAMPLE:
     559        sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import all_labeled_graphs
     560        sage: st = sage.groups.perm_gps.partn_ref.refinement_graphs.search_tree
     561        sage: Glist = {}
     562        sage: Giso  = {}
     563        sage: for n in [1..5]:
     564        ...    Glist[n] = all_labeled_graphs(n)
     565        ...    Giso[n] = []
     566        ...    for g in Glist[n]:
     567        ...        a, b = st(g, [range(n)])
     568        ...        inn = False
     569        ...        for gi in Giso[n]:
     570        ...            if b == gi:
     571        ...                inn = True
     572        ...        if not inn:
     573        ...            Giso[n].append(b)
     574        sage: for n in Giso:
     575        ...    print n, len(Giso[n])
     576        1 1
     577        2 2
     578        3 4
     579        4 11
     580        5 34
     581
     582    """
     583    from sage.graphs.graph import Graph
     584    TE = []
     585    for i in range(n):
     586        for j in range(i):
     587            TE.append((i, j))
     588    m = len(TE)
     589    Glist= []
     590    for i in range(2**m):
     591        G = Graph(n)
     592        b = Integer(i).binary()
     593        b = '0'*(m-len(b)) + b
     594        for i in range(m):
     595            if int(b[i]):
     596                G.add_edge(TE[i])
     597        Glist.append(G)
     598    return Glist
     599
     600
     601def random_tests(t=10.0, n_max=60, perms_per_graph=10):
     602    """
     603    Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma
     604    and random graphs G.
     605   
     606    DOCTEST:
     607        sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
     608        sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests()
     609        All passed: ... random tests on ... graphs.
     610        sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests(180.0, 200, 30) # long time
     611        All passed: ... random tests on ... graphs.       
     612
     613    """
     614    from sage.misc.misc import walltime
     615    from sage.misc.prandom import random, randint
     616    from sage.graphs.graph_generators import GraphGenerators, DiGraphGenerators
     617    from sage.combinat.permutation import Permutations
     618    cdef int i, j, num_tests = 0, num_graphs = 0, passed = 1
     619    GG = GraphGenerators()
     620    DGG = DiGraphGenerators()
     621    t_0 = walltime()
     622    while walltime(t_0) < t:
     623        p = random()
     624        n = randint(1, n_max)
     625        S = Permutations(n)
     626
     627        G = GG.RandomGNP(n, p)
     628        H = G.copy()
     629        for i from 0 <= i < perms_per_graph:
     630            G1 = search_tree(G, [G.vertices()])[1]
     631            perm = list(S.random_element())
     632            perm = [perm[j]-1 for j from 0 <= j < n]
     633            G.relabel(perm)
     634            G2 = search_tree(G, [G.vertices()])[1]
     635            if G1 != G2:
     636                print "FAILURE: graph6-"
     637                print H.graph6_string()
     638                print perm
     639                passed = 0
     640
     641        D = DGG.RandomDirectedGNP(n, p)
     642        D.loops(True)
     643        for i from 0 <= i < n:
     644            if random() < p:
     645                D.add_edge(i,i)
     646        E = D.copy()
     647        for i from 0 <= i < perms_per_graph:
     648            D1 = search_tree(D, [D.vertices()], dig=True)[1]
     649            perm = list(S.random_element())
     650            perm = [perm[j]-1 for j from 0 <= j < n]
     651            D.relabel(perm)
     652            D2 = search_tree(D, [D.vertices()], dig=True)[1]
     653            if D1 != D2:
     654                print "FAILURE: dig6-"
     655                print E.dig6_string()
     656                print perm
     657                passed = 0
     658        num_tests += 2*perms_per_graph
     659        num_graphs += 2
     660    if passed:
     661        print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs)
     662
     663
  • new file sage/groups/perm_gps/partn_ref/tree_traversal.pxd

    diff -r 91af4c3f6b92 -r 074731a26986 sage/groups/perm_gps/partn_ref/tree_traversal.pxd
    - +  
     1
     2#*****************************************************************************
     3#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     4#
     5# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     6#                         http://www.gnu.org/licenses/
     7#*****************************************************************************
     8
     9include '../../../ext/cdefs.pxi'
     10include '../../../ext/stdsage.pxi'
     11from sage.rings.integer cimport Integer
     12
     13cdef struct OrbitPartition:
     14    # Disjoint set data structure for representing the orbits of the generators
     15    # so far found. Also keeps track of the minimum elements of the cells and
     16    # their sizes.
     17    int degree
     18    int *parent
     19    int *rank
     20    int *mcr # minimum cell representatives - only valid at the root of a cell
     21    int *size # also only valid at the root of a cell
     22
     23cdef inline OrbitPartition *OP_new(int)
     24cdef inline int OP_dealloc(OrbitPartition *)
     25cdef inline int OP_find(OrbitPartition *, int)
     26cdef inline int OP_join(OrbitPartition *, int, int)
     27cdef inline int OP_merge_list_perm(OrbitPartition *, int *)
     28
     29cdef struct PartitionStack:
     30    # Representation of a node of the search tree. A sequence of partitions of
     31    # length depth + 1, each of which is finer than the last. Partition k is
     32    # represented as PS.entries in order, broken up immediately after each
     33    # entry of levels which is at most k.
     34    int *entries
     35    int *levels
     36    int depth
     37    int degree
     38
     39cdef inline PartitionStack *PS_new(int, bint)
     40cdef inline PartitionStack *PS_copy(PartitionStack *)
     41cdef inline int PS_copy_from_to(PartitionStack *, PartitionStack *)
     42cdef inline PartitionStack *PS_from_list(object)
     43cdef inline int PS_dealloc(PartitionStack *)
     44cdef inline int PS_print(PartitionStack *)
     45cdef inline int PS_print_partition(PartitionStack *, int)
     46cdef inline bint PS_is_discrete(PartitionStack *)
     47cdef inline int PS_num_cells(PartitionStack *)
     48cdef inline int PS_move_min_to_front(PartitionStack *, int, int)
     49cdef inline bint PS_is_mcr(PartitionStack *, int)
     50cdef inline bint PS_is_fixed(PartitionStack *, int)
     51cdef inline int PS_clear(PartitionStack *)
     52cdef inline int PS_split_point(PartitionStack *, int)
     53#cdef inline int PS_first_smallest(PartitionStack *, bitset_t)
     54cdef inline int PS_get_perm_from(PartitionStack *, PartitionStack *, int *)
     55
     56cdef inline int split_point_and_refine(PartitionStack *, int, object,
     57    int (*)(PartitionStack *, object, int *, int), int *)
     58
     59cdef struct traverse_return:
     60    int *generators
     61    int num_gens
     62    int *relabeling
     63    int *base
     64    int base_size
     65    mpz_t order
     66
     67cdef traverse_return *traverse_tree( object, int **, int,
     68    bint (*)(PartitionStack *, object),
     69    int (*)(PartitionStack *, object, int *, int),
     70    int (*)(int *, int *, object), bint, bint, bint)
     71
     72
     73
  • new file sage/groups/perm_gps/partn_ref/tree_traversal.pyx

    diff -r 91af4c3f6b92 -r 074731a26986 sage/groups/perm_gps/partn_ref/tree_traversal.pyx
    - +  
     1"""
     2Module for traversing the search tree associated with partition backtrack
     3problems.
     4
     5DOCTEST:
     6    sage: import sage.groups.perm_gps.partn_ref.tree_traversal
     7
     8"""
     9
     10#*****************************************************************************
     11#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     12#
     13# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     14#                         http://www.gnu.org/licenses/
     15#*****************************************************************************
     16
     17include '../../../misc/bitset.pxi'
     18
     19cdef inline int min(int a, int b):
     20    if a < b: return a
     21    else: return b
     22
     23cdef inline int max(int a, int b):
     24    if a > b: return a
     25    else: return b
     26
     27cdef inline int cmp(int a, int b):
     28    if a < b: return -1
     29    elif a == b: return 0
     30    else: return 1
     31
     32# OrbitPartitions
     33
     34cdef inline OrbitPartition *OP_new(int n):
     35    """
     36    Allocate and return a pointer to a new OrbitPartition of degree n. Returns a
     37    null pointer in the case of an allocation failure.
     38    """
     39    cdef int i
     40    cdef OrbitPartition *OP = <OrbitPartition *> \
     41                                sage_malloc(sizeof(OrbitPartition))
     42    if OP is NULL:
     43        return OP
     44    OP.degree = n
     45    OP.parent = <int *> sage_malloc( n * sizeof(int) )
     46    OP.rank = <int *> sage_malloc( n * sizeof(int) )
     47    OP.mcr = <int *> sage_malloc( n * sizeof(int) )
     48    OP.size = <int *> sage_malloc( n * sizeof(int) )
     49    if OP.parent is NULL or OP.rank is NULL or \
     50       OP.mcr is NULL or OP.size is NULL:
     51        if OP.parent is not NULL: sage_free(OP.parent)
     52        if OP.rank is not NULL: sage_free(OP.rank)
     53        if OP.mcr is not NULL: sage_free(OP.mcr)
     54        if OP.size is not NULL: sage_free(OP.size)
     55        return NULL
     56    for i from 0 <= i < n:
     57        OP.parent[i] = i
     58        OP.rank[i] = 0
     59        OP.mcr[i] = i
     60        OP.size[i] = 1
     61    return OP
     62
     63cdef inline int OP_dealloc(OrbitPartition *OP):
     64    sage_free(OP.parent)
     65    sage_free(OP.rank)
     66    sage_free(OP.mcr)
     67    sage_free(OP.size)
     68    sage_free(OP)
     69    return 0
     70
     71cdef inline int OP_find(OrbitPartition *OP, int n):
     72    """
     73    Report the representative ("root") of the cell which contains n.
     74    """
     75    if OP.parent[n] == n:
     76        return n
     77    else:
     78        OP.parent[n] = OP_find(OP, OP.parent[n])
     79        return OP.parent[n]
     80
     81cdef inline int OP_join(OrbitPartition *OP, int m, int n):
     82    """
     83    Join the cells containing m and n, if they are different.
     84    """
     85    cdef int m_root = OP_find(OP, m)
     86    cdef int n_root = OP_find(OP, n)
     87    if OP.rank[m_root] > OP.rank[n_root]:
     88        OP.parent[n_root] = m_root
     89        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
     90        OP.size[m_root] += OP.size[n_root]
     91    elif OP.rank[m_root] < OP.rank[n_root]:
     92        OP.parent[m_root] = n_root
     93        OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root])
     94        OP.size[n_root] += OP.size[m_root]
     95    elif m_root != n_root:
     96        OP.parent[n_root] = m_root
     97        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
     98        OP.size[m_root] += OP.size[n_root]
     99        OP.rank[m_root] += 1
     100
     101cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma):
     102    """
     103    Joins the cells of OP which intersect the same orbit of gamma.
     104   
     105    INPUT:
     106        gamma - an integer array representing i -> gamma[i].
     107   
     108    OUTPUT:
     109        1 - something changed
     110        0 - orbits of gamma all contained in cells of OP
     111    """
     112    cdef int i, i_root, gamma_i_root, changed = 0
     113    for i from 0 <= i < OP.degree:
     114        if gamma[i] == i: continue
     115        i_root = OP_find(OP, i)
     116        gamma_i_root = OP_find(OP, gamma[i])
     117        if i_root != gamma_i_root:
     118            changed = 1
     119            OP_join(OP, i_root, gamma_i_root)
     120    return changed
     121
     122def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]):
     123    """
     124    Demonstration and testing.
     125   
     126    DOCTEST:
     127        sage: from sage.groups.perm_gps.partn_ref.tree_traversal \
     128        ... import OP_represent
     129        sage: OP_represent()
     130        Allocating OrbitPartition...
     131        Allocation passed.
     132        Checking that each element reports itself as its root.
     133        Each element reports itself as its root.
     134        Merging:
     135        Merged 0 and 1.
     136        Merged 2 and 3.
     137        Merged 3 and 4.
     138        Done merging.
     139        Finding:
     140        0 -> 0, root: size=2, mcr=0, rank=1
     141        1 -> 0
     142        2 -> 2, root: size=3, mcr=2, rank=1
     143        3 -> 2
     144        4 -> 2
     145        5 -> 5, root: size=1, mcr=5, rank=0
     146        6 -> 6, root: size=1, mcr=6, rank=0
     147        7 -> 7, root: size=1, mcr=7, rank=0
     148        8 -> 8, root: size=1, mcr=8, rank=0
     149        Allocating array to test merge_perm.
     150        Allocation passed.
     151        Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8]
     152        Done merging.
     153        Finding:
     154        0 -> 0, root: size=5, mcr=0, rank=2
     155        1 -> 0
     156        2 -> 0
     157        3 -> 0
     158        4 -> 0
     159        5 -> 5, root: size=3, mcr=5, rank=1
     160        6 -> 5
     161        7 -> 5
     162        8 -> 8, root: size=1, mcr=8, rank=0
     163        Deallocating OrbitPartition.
     164        Done.
     165
     166    """
     167    cdef int i
     168    print "Allocating OrbitPartition..."
     169    cdef OrbitPartition *OP = OP_new(n)
     170    if OP is NULL:
     171        print "Allocation failed!"
     172        return
     173    print "Allocation passed."
     174    print "Checking that each element reports itself as its root."
     175    good = True
     176    for i from 0 <= i < n:
     177        if not OP_find(OP, i) == i:
     178            print "Failed at i = %d!"%i
     179            good = False
     180    if good: print "Each element reports itself as its root."
     181    print "Merging:"
     182    for i,j in merges:
     183        OP_join(OP, i, j)
     184        print "Merged %d and %d."%(i,j)
     185    print "Done merging."
     186    print "Finding:"
     187    for i from 0 <= i < n:
     188        j = OP_find(OP, i)
     189        s = "%d -> %d"%(i, j)
     190        if i == j:
     191            s += ", root: size=%d, mcr=%d, rank=%d"%\
     192                   (OP.size[i], OP.mcr[i], OP.rank[i])
     193        print s
     194    print "Allocating array to test merge_perm."
     195    cdef int *gamma = <int *> sage_malloc( n * sizeof(int) )
     196    if gamma is NULL:
     197        print "Allocation failed!"
     198        OP_dealloc(OP)
     199        return
     200    print "Allocation passed."
     201    for i from 0 <= i < n:
     202        gamma[i] = perm[i]
     203    print "Merging permutation: %s"%perm
     204    OP_merge_list_perm(OP, gamma)
     205    print "Done merging."
     206    print "Finding:"
     207    for i from 0 <= i < n:
     208        j = OP_find(OP, i)
     209        s = "%d -> %d"%(i, j)
     210        if i == j:
     211            s += ", root: size=%d, mcr=%d, rank=%d"%\
     212                   (OP.size[i], OP.mcr[i], OP.rank[i])
     213        print s   
     214    print "Deallocating OrbitPartition."
     215    sage_free(gamma)
     216    OP_dealloc(OP)
     217    print "Done."
     218
     219# PartitionStacks
     220
     221cdef inline PartitionStack *PS_new(int n, bint unit_partition):
     222    """
     223    Allocate and return a pointer to a new PartitionStack of degree n. Returns a
     224    null pointer in the case of an allocation failure.
     225    """
     226    cdef int i
     227    cdef PartitionStack *PS = <PartitionStack *> \
     228                                sage_malloc(sizeof(PartitionStack))
     229    if PS is NULL:
     230        return PS
     231    PS.entries = <int *> sage_malloc( n * sizeof(int) )
     232    PS.levels = <int *> sage_malloc( n * sizeof(int) )
     233    if PS.entries is NULL or PS.levels is NULL:
     234        if PS.entries is not NULL: sage_free(PS.entries)
     235        if PS.levels is not NULL: sage_free(PS.levels)
     236        return NULL
     237    PS.depth = 0
     238    PS.degree = n
     239    if unit_partition:
     240        for i from 0 <= i < n-1:
     241            PS.entries[i] = i
     242            PS.levels[i] = n
     243        PS.entries[n-1] = n-1
     244        PS.levels[n-1] = -1
     245    return PS
     246
     247cdef inline PartitionStack *PS_copy(PartitionStack *PS):
     248    """
     249    Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
     250    pointer in the case of an allocation failure.
     251    """
     252    cdef int i
     253    cdef PartitionStack *PS2 = <PartitionStack *> \
     254                                sage_malloc(sizeof(PartitionStack))
     255    if PS2 is NULL:
     256        return PS2
     257    PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
     258    PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
     259    if PS2.entries is NULL or PS2.levels is NULL:
     260        if PS2.entries is not NULL: sage_free(PS2.entries)
     261        if PS2.levels is not NULL: sage_free(PS2.levels)
     262        return NULL
     263    PS_copy_from_to(PS, PS2)
     264    return PS2
     265
     266cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2):
     267    """
     268    Copy all data from PS to PS2.
     269    """
     270    PS2.depth = PS.depth
     271    PS2.degree = PS.degree
     272    for i from 0 <= i < PS.degree:
     273        PS2.entries[i] = PS.entries[i]
     274        PS2.levels[i] = PS.levels[i]   
     275
     276cdef inline PartitionStack *PS_from_list(object L):
     277    """
     278    Allocate and return a pointer to a PartitionStack representing L. Returns a
     279    null pointer in the case of an allocation failure.
     280    """
     281    cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0
     282    for cell from 0 <= cell < num_cells:
     283        n += len(L[cell])
     284    cdef PartitionStack *PS = PS_new(n, False)
     285    cell = 0
     286    if PS is NULL:
     287        return PS
     288    while 1:
     289        cur_len = len(L[cell])
     290        for i from 0 <= i < cur_len:
     291            PS.entries[cur_start + i] = L[cell][i]
     292            PS.levels[cur_start + i] = n
     293        PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
     294        cur_start += cur_len
     295        cell += 1
     296        if cell == num_cells:
     297            PS.levels[cur_start-1] = -1
     298            break
     299        PS.levels[cur_start-1] = 0
     300    return PS
     301
     302cdef inline int PS_dealloc(PartitionStack *PS):
     303    sage_free(PS.entries)
     304    sage_free(PS.levels)
     305    sage_free(PS)
     306    return 0
     307
     308cdef inline int PS_print(PartitionStack *PS):
     309    """
     310    Print a visual representation of PS.
     311    """
     312    cdef int i
     313    for i from 0 <= i < PS.depth + 1:
     314        PS_print_partition(PS, i)
     315
     316cdef inline int PS_print_partition(PartitionStack *PS, int k):
     317    """
     318    Print the partition at depth k.
     319    """
     320    s = '('
     321    for i from 0 <= i < PS.degree:
     322        s += str(PS.entries[i])
     323        if PS.levels[i] <= k:
     324            s += '|'
     325        else:
     326            s += ','
     327    s = s[:-1] + ')'
     328    print s
     329
     330cdef inline bint PS_is_discrete(PartitionStack *PS):
     331    """
     332    Returns whether the deepest partition consists only of singleton cells.
     333    """
     334    cdef int i
     335    for i from 0 <= i < PS.degree:
     336        if PS.levels[i] > PS.depth:
     337            return 0
     338    return 1
     339
     340cdef inline int PS_num_cells(PartitionStack *PS):
     341    """
     342    Returns the number of cells.
     343    """
     344    cdef int i, ncells = 0
     345    for i from 0 <= i < PS.degree:
     346        if PS.levels[i] <= PS.depth:
     347            ncells += 1
     348    return ncells
     349
     350cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
     351    """
     352    Makes sure that the first element of the segment of entries i with
     353    start <= i <= end is minimal.
     354    """
     355    cdef int i, min_loc = start, minimum = PS.entries[start]
     356    for i from start < i <= end:
     357        if PS.entries[i] < minimum:
     358            min_loc = i
     359            minimum = PS.entries[i]
     360    if min_loc != start:
     361        PS.entries[min_loc] = PS.entries[start]
     362        PS.entries[start] = minimum
     363
     364cdef inline bint PS_is_mcr(PartitionStack *PS, int m):
     365    """
     366    Returns whether PS.elements[m] (not m!) is the smallest element of its cell.
     367    """
     368    return m == 0 or PS.levels[m-1] <= PS.depth
     369
     370cdef inline bint PS_is_fixed(PartitionStack *PS, int m):
     371    """
     372    Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming
     373    PS_is_mcr(PS, m) is already true.
     374    """
     375    return PS.levels[m] <= PS.depth
     376
     377cdef inline int PS_clear(PartitionStack *PS):
     378    """
     379    Sets the current partition to the first shallower one, i.e. forgets about
     380    boundaries between cells that are new to the current level.
     381    """
     382    cdef int i, cur_start = 0
     383    for i from 0 <= i < PS.degree:
     384        if PS.levels[i] >= PS.depth:
     385            PS.levels[i] += 1
     386        if PS.levels[i] < PS.depth:
     387            PS_move_min_to_front(PS, cur_start, i)
     388            cur_start = i+1
     389
     390cdef inline int PS_split_point(PartitionStack *PS, int v):
     391    """
     392    Detaches the point v from the cell it is in, putting the singleton cell
     393    of just v in front. Returns the position where v is now located.
     394    """
     395    cdef int i = 0, index_of_v
     396    while PS.entries[i] != v:
     397        i += 1
     398    index_of_v = i
     399    while PS.levels[i] > PS.depth:
     400        i += 1
     401    if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \
     402       and PS.levels[index_of_v] > PS.depth:
     403        # if v is first (make sure v is not already alone)
     404        PS_move_min_to_front(PS, index_of_v+1, i)
     405        PS.levels[index_of_v] = PS.depth
     406        return index_of_v
     407    else:
     408        # If v is not at front, i.e. v is not minimal in its cell,
     409        # then move_min_to_front is not necessary since v will swap
     410        # with the first before making its own cell, leaving it at
     411        # the front of the other.
     412        i = index_of_v
     413        while i != 0 and PS.levels[i-1] > PS.depth:
     414            i -= 1
     415        PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v
     416        PS.entries[i+1] = PS.entries[i] # move the first (min) to second
     417        PS.entries[i] = v # place v first
     418        PS.levels[i] = PS.depth
     419        return i
     420
     421cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
     422    """
     423    Find the first occurrence of the smallest cell of size greater than one,
     424    store its entries to b, and return its minimum element.
     425    """
     426    cdef int i = 0, j = 0, location = 0, n = PS.degree
     427    bitset_zero(b)
     428    while 1:
     429        if PS.levels[i] <= PS.depth:
     430            if i != j and n > i - j + 1:
     431                n = i - j + 1
     432                location = j
     433            j = i + 1
     434        if PS.levels[i] == -1: break
     435        i += 1
     436    # location now points to the beginning of the first, smallest,
     437    # nontrivial cell
     438    i = location
     439    while 1:
     440        bitset_flip(b, PS.entries[i])
     441        if PS.levels[i] <= PS.depth: break
     442        i += 1
     443    return PS.entries[location]
     444
     445cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
     446    """
     447    Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
     448    denotes the entry of the ith cell of the discrete partition PS.
     449    """
     450    cdef int i = 0
     451    for i from 0 <= i < PS1.degree:
     452        gamma[PS2.entries[i]] = PS1.entries[i]
     453
     454def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]):
     455    """
     456    Demonstration and testing.
     457   
     458    DOCTEST:
     459        sage: from sage.groups.perm_gps.partn_ref.tree_traversal \
     460        ... import PS_represent
     461        sage: PS_represent()
     462        Allocating PartitionStack...
     463        Allocation passed:
     464        (0,1,2,3,4,5,6,7,8,9)
     465        Checking that entries are in order and correct level.
     466        Everything seems in order, deallocating.
     467        Deallocated.
     468        Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]].
     469        PartitionStack's data:
     470        entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2]
     471        levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1]
     472        depth = 0, degree = 10
     473        (6|3,4,8,7|1,9,5|0,2)
     474        Checking PS_is_discrete:
     475        False
     476        Checking PS_num_cells:
     477        4
     478        Checking PS_is_mcr, min cell reps are:
     479        [6, 3, 1, 0]
     480        Checking PS_is_fixed, fixed elements are:
     481        [6]
     482        Copying PartitionStack:
     483        (6|3,4,8,7|1,9,5|0,2)
     484        Checking for consistency.
     485        Everything is consistent.
     486        Clearing copy:
     487        (0,3,4,8,7,1,9,5,6,2)
     488        Splitting point 6 from original:
     489        0
     490        (6|3,4,8,7|1,9,5|0,2)
     491        Splitting point 1 from original:
     492        5
     493        (6|3,4,8,7|1|5,9|0,2)
     494        Splitting point 8 from original:
     495        1
     496        (6|8|3,4,7|1|5,9|0,2)
     497        Splitting point 2 from original:
     498        8
     499        (6|8|3,4,7|1|5,9|2|0)
     500        Getting permutation from PS2->PS:
     501        [6, 1, 0, 8, 3, 9, 2, 7, 4, 5]
     502        Finding first smallest:
     503        Minimal element is 5, bitset is:
     504        0000010001
     505        Deallocating PartitionStacks.
     506        Done.
     507
     508    """
     509    cdef int i, n = sum([len(cell) for cell in partition]), *gamma
     510    cdef bitset_t b
     511    print "Allocating PartitionStack..."
     512    cdef PartitionStack *PS = PS_new(n, True), *PS2
     513    if PS is NULL:
     514        print "Allocation failed!"
     515        return
     516    print "Allocation passed:"
     517    PS_print(PS)
     518    print "Checking that entries are in order and correct level."
     519    good = True
     520    for i from 0 <= i < n-1:
     521        if not (PS.entries[i] == i and PS.levels[i] == n):
     522            print "Failed at i = %d!"%i
     523            print PS.entries[i], PS.levels[i], i, n
     524            good = False
     525    if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1):
     526        print "Failed at i = %d!"%(n-1)
     527        good = False
     528    if not PS.degree == n or not PS.depth == 0:
     529        print "Incorrect degree or depth!"
     530        good = False
     531    if good: print "Everything seems in order, deallocating."
     532    PS_dealloc(PS)
     533    print "Deallocated."
     534    print "Creating PartitionStack from partition %s."%partition
     535    PS = PS_from_list(partition)
     536    print "PartitionStack's data:"
     537    print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n]
     538    print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n]
     539    print "depth = %d, degree = %d"%(PS.depth,PS.degree)
     540    PS_print(PS)
     541    print "Checking PS_is_discrete:"
     542    print "True" if PS_is_discrete(PS) else "False"
     543    print "Checking PS_num_cells:"
     544    print PS_num_cells(PS)
     545    print "Checking PS_is_mcr, min cell reps are:"
     546    L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)]
     547    print L
     548    print "Checking PS_is_fixed, fixed elements are:"
     549    print [PS.entries[l] for l in L if PS_is_fixed(PS, l)]
     550    print "Copying PartitionStack:"
     551    PS2 = PS_copy(PS)
     552    PS_print(PS2)
     553    print "Checking for consistency."
     554    good = True
     555    for i from 0 <= i < n:
     556        if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]:
     557            print "Failed at i = %d!"%i
     558            good = False
     559    if PS.degree != PS2.degree or PS.depth != PS2.depth:
     560        print "Failure with degree or depth!"
     561        good = False
     562    if good:
     563        print "Everything is consistent."
     564    print "Clearing copy:"
     565    PS_clear(PS2)
     566    PS_print(PS2)
     567    for s in splits:
     568        print "Splitting point %d from original:"%s
     569        print PS_split_point(PS, s)
     570        PS_print(PS)
     571    print "Getting permutation from PS2->PS:"
     572    gamma = <int *> sage_malloc(n * sizeof(int))
     573    PS_get_perm_from(PS, PS2, gamma)
     574    print [gamma[i] for i from 0 <= i < n]
     575    sage_free(gamma)
     576    print "Finding first smallest:"
     577    bitset_init(b, n)
     578    i = PS_first_smallest(PS, b)
     579    print "Minimal element is %d, bitset is:"%i
     580    print bitset_string(b)
     581    bitset_clear(b)
     582    print "Deallocating PartitionStacks."
     583    PS_dealloc(PS)
     584    PS_dealloc(PS2)
     585    print "Done."
     586
     587# Functions
     588
     589cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
     590    int (*refine_and_return_invariant)\
     591         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
     592    int *cells_to_refine_by):
     593    """
     594    Make the partition stack one longer by copying the last partition in the
     595    stack, split off a given point, and refine. Return the invariant given by
     596    the refinement function.
     597   
     598    INPUT:
     599    PS -- the partition stack to refine
     600    v -- the point to split
     601    S -- the structure
     602    refine_and_return_invariant -- the refinement function provided
     603    cells_to_refine_by -- an array, contents ignored
     604
     605    """
     606    PS.depth += 1
     607    PS_clear(PS)
     608    cells_to_refine_by[0] = PS_split_point(PS, v)
     609    return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
     610
     611cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, object S):
     612    return 0
     613
     614cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
     615    return 0
     616
     617cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S):
     618    return 0
     619
     620def test_traverse_tree_trivially(int n=6, partition=[[0,1,2],[3,4],[5]],
     621    canonical_label=True, base=False):
     622    """
     623    sage: tttt = sage.groups.perm_gps.partn_ref.tree_traversal.test_traverse_tree_trivially
     624    sage: tttt()
     625    12
     626    sage: tttt(canonical_label=False, base=False)
     627    12
     628    sage: tttt(canonical_label=False, base=True)
     629    12
     630    sage: tttt(canonical_label=True, base=True)
     631    12
     632    sage: tttt(n=0, partition=[])
     633    1
     634    sage: tttt(n=0, partition=[], canonical_label=False, base=False)
     635    1
     636    sage: tttt(n=0, partition=[], canonical_label=False, base=True)
     637    1
     638    sage: tttt(n=0, partition=[], canonical_label=True, base=True)
     639    1
     640   
     641    """
     642    cdef int i, j
     643    cdef traverse_return *output
     644    cdef Integer I = Integer(0)
     645    cdef int **part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     646    for i from 0 <= i < len(partition):
     647        part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
     648        for j from 0 <= j < len(partition[i]):
     649            part[i][j] = partition[i][j]
     650        part[i][len(partition[i])] = -1
     651    part[len(partition)] = NULL
     652    output = traverse_tree([], part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_structures_trivial, canonical_label, base, True)
     653    mpz_set(I.value, output.order)
     654    print I
     655    for i from 0 <= i < len(partition):
     656        sage_free(part[i])
     657    sage_free(part)
     658    mpz_clear(output.order)
     659    sage_free(output.generators)
     660    if base:
     661        sage_free(output.base)
     662    if canonical_label:
     663        sage_free(output.relabeling)
     664    sage_free(output)
     665
     666cdef traverse_return *traverse_tree(object S, int **partition, int n,
     667    bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
     668    int (*refine_and_return_invariant)\
     669         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
     670    int (*compare_structures)(int *gamma_1, int *gamma_2, object S),
     671    bint canonical_label, bint base, bint order):
     672    """
     673    Traverse the search space for subgroup/canonical label calculation.
     674   
     675    INPUT:
     676    S -- pointer to the structure
     677    partition -- array representing a partition of the points
     678    len_partition -- length of the partition
     679    n -- the number of points (points are assumed to be 0,1,...,n-1)
     680    canonical_label -- whether to search for canonical label - if True, return
     681        the permutation taking S to its canonical label
     682    base -- if True, return the base for which the given generating set is a
     683        strong generating set
     684    order -- if True, return the order of the automorphism group
     685    all_children_are_equivalent -- pointer to a function
     686        INPUT:
     687        PS -- pointer to a partition stack
     688        S -- pointer to the structure
     689        OUTPUT:
     690        bint -- returns True if it can be determined that all refinements below
     691            the current one will result in an equivalent discrete partition
     692    refine_and_return_invariant -- pointer to a function
     693        INPUT:
     694        PS -- pointer to a partition stack
     695        S -- pointer to the structure
     696        alpha -- an array consisting of numbers, which indicate the starting
     697            positions of the cells to refine against (will likely be modified)
     698        OUTPUT:
     699        int -- returns an invariant under application of arbitrary permutations
     700    compare_structures -- pointer to a function
     701        INPUT:
     702        gamma_1, gamma_2 -- (list) permutations of the points of S
     703        S -- pointer to the structure
     704        OUTPUT:
     705        int -- 0 if gamma_1(S) = gamma_2(S), otherwise -1 or 1 (see docs for cmp),
     706            such that the set of all structures is well-ordered
     707    OUTPUT:
     708    pointer to a traverse_return struct
     709
     710    """
     711    cdef PartitionStack *current_ps, *first_ps, *label_ps
     712    cdef int first_meets_current = -1
     713    cdef int label_meets_current
     714    cdef int current_kids_are_same = 1
     715    cdef int first_kids_are_same
     716   
     717    cdef int *current_indicators, *first_indicators, *label_indicators
     718    cdef int first_and_current_indicator_same
     719    cdef int label_and_current_indicator_same = -1
     720    cdef int compared_current_and_label_indicators
     721   
     722    cdef OrbitPartition *orbits_of_subgroup
     723    cdef mpz_t subgroup_size
     724    cdef int subgroup_primary_orbit_size = 0
     725    cdef int minimal_in_primary_orbit
     726    cdef int size_of_generator_array = 2*n*n
     727   
     728    cdef bitset_t *fixed_points_of_generators # i.e. fp
     729    cdef bitset_t *minimal_cell_reps_of_generators # i.e. mcr
     730    cdef int len_of_fp_and_mcr = 100
     731    cdef int index_in_fp_and_mcr = -1
     732   
     733    cdef bitset_t *vertices_to_split
     734    cdef bitset_t vertices_have_been_reduced
     735    cdef int *permutation, *id_perm, *cells_to_refine_by
     736    cdef int *vertices_determining_current_stack
     737   
     738    cdef int i, j, k
     739    cdef bint discrete, automorphism, update_label
     740    cdef bint backtrack, new_vertex, narrow, mem_err = 0
     741   
     742    cdef traverse_return *output = <traverse_return *> sage_malloc(sizeof(traverse_return))
     743    if output is NULL:
     744        raise MemoryError
     745   
     746    if n == 0:
     747        output.generators = NULL
     748        output.num_gens = 0
     749        output.relabeling = NULL
     750        output.base = NULL
     751        output.base_size = 0
     752        mpz_init_set_si(output.order, 1)
     753        return output
     754   
     755    # Allocate:
     756    mpz_init_set_si(subgroup_size, 1)
     757   
     758    output.generators = <int *> sage_malloc( size_of_generator_array * sizeof(int) )
     759    output.num_gens = 0
     760    if base:
     761        output.base = <int *> sage_malloc( n * sizeof(int) )
     762   
     763    current_indicators = <int *> sage_malloc( n * sizeof(int) )
     764    first_indicators = <int *> sage_malloc( n * sizeof(int) )
     765   
     766    if canonical_label:
     767        label_indicators = <int *> sage_malloc( n * sizeof(int) )
     768        output.relabeling = <int *> sage_malloc( n * sizeof(int) )
     769
     770    fixed_points_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
     771    minimal_cell_reps_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
     772   
     773    vertices_to_split = <bitset_t *> sage_malloc( n * sizeof(bitset_t) )
     774    permutation = <int *> sage_malloc( n * sizeof(int) )
     775    id_perm = <int *> sage_malloc( n * sizeof(int) )
     776    cells_to_refine_by = <int *> sage_malloc( n * sizeof(int) )
     777    vertices_determining_current_stack = <int *> sage_malloc( n * sizeof(int) )
     778   
     779    current_ps = PS_new(n, False)
     780    orbits_of_subgroup = OP_new(n)
     781   
     782    # Check for allocation failures:
     783    cdef bint succeeded = 1
     784    if canonical_label:
     785        if label_indicators is NULL or output.relabeling is NULL:
     786            succeeded = 0
     787            if label_indicators is not NULL:
     788                sage_free(label_indicators)
     789            if output.relabeling is not NULL:
     790                sage_free(output.relabeling)
     791    if base:
     792        if output.base is NULL:
     793            succeeded = 0
     794    if not succeeded                              or \
     795       output.generators                  is NULL or \
     796       current_indicators                 is NULL or \
     797       first_indicators                   is NULL or \
     798       fixed_points_of_generators         is NULL or \
     799       minimal_cell_reps_of_generators    is NULL or \
     800       vertices_to_split                  is NULL or \
     801       permutation                        is NULL or \
     802       id_perm                            is NULL or \
     803       cells_to_refine_by                 is NULL or \
     804       vertices_determining_current_stack is NULL or \
     805       current_ps                         is NULL or \
     806       orbits_of_subgroup                 is NULL:
     807        if output.generators                  is not NULL:
     808            sage_free(output.generators)
     809        if base:
     810            if output.base                        is not NULL:
     811                sage_free(output.base)
     812        if current_indicators                 is not NULL:
     813            sage_free(current_indicators)
     814        if first_indicators                   is not NULL:
     815            sage_free(first_indicators)
     816        if fixed_points_of_generators         is not NULL:
     817            sage_free(fixed_points_of_generators)
     818        if minimal_cell_reps_of_generators    is not NULL:
     819            sage_free(minimal_cell_reps_of_generators)
     820        if vertices_to_split                  is not NULL:
     821            sage_free(vertices_to_split)
     822        if permutation                        is not NULL:
     823            sage_free(permutation)
     824        if id_perm                            is not NULL:
     825            sage_free(id_perm)
     826        if cells_to_refine_by                 is not NULL:
     827            sage_free(cells_to_refine_by)
     828        if vertices_determining_current_stack is not NULL:
     829            sage_free(vertices_determining_current_stack)
     830        if current_ps is not NULL:
     831            PS_dealloc(current_ps)
     832        if orbits_of_subgroup is not NULL:
     833            OP_dealloc(orbits_of_subgroup)
     834        sage_free(output)
     835        mpz_clear(subgroup_size)
     836        raise MemoryError
     837   
     838    # Initialize bitsets, checking for allocation failures:
     839    succeeded = 1
     840    for i from 0 <= i < len_of_fp_and_mcr:
     841        try:
     842            bitset_init(fixed_points_of_generators[i], n)
     843        except MemoryError:
     844            succeeded = 0
     845            for j from 0 <= j < i:
     846                bitset_clear(fixed_points_of_generators[j])
     847                bitset_clear(minimal_cell_reps_of_generators[j])
     848            break
     849        try:
     850            bitset_init(minimal_cell_reps_of_generators[i], n)
     851        except MemoryError:
     852            succeeded = 0
     853            for j from 0 <= j < i:
     854                bitset_clear(fixed_points_of_generators[j])
     855                bitset_clear(minimal_cell_reps_of_generators[j])
     856            bitset_clear(fixed_points_of_generators[i])
     857            break
     858    if succeeded:
     859        for i from 0 <= i < n:
     860            try:
     861                bitset_init(vertices_to_split[i], n)
     862            except MemoryError:
     863                succeeded = 0
     864                for j from 0 <= j < i:
     865                    bitset_clear(vertices_to_split[j])
     866                for j from 0 <= j < len_of_fp_and_mcr:
     867                    bitset_clear(fixed_points_of_generators[j])
     868                    bitset_clear(minimal_cell_reps_of_generators[j])
     869                break
     870    if succeeded:
     871        try:
     872            bitset_init(vertices_have_been_reduced, n)
     873        except MemoryError:
     874            succeeded = 0
     875            for j from 0 <= j < n:
     876                bitset_clear(vertices_to_split[j])
     877            for j from 0 <= j < len_of_fp_and_mcr:
     878                bitset_clear(fixed_points_of_generators[j])
     879                bitset_clear(minimal_cell_reps_of_generators[j])
     880    if not succeeded:
     881        sage_free(output.generators)
     882        if base:
     883            sage_free(output.base)
     884        sage_free(current_indicators)
     885        sage_free(first_indicators)
     886        if canonical_label:
     887            sage_free(output.relabeling)
     888            sage_free(label_indicators)
     889        sage_free(fixed_points_of_generators)
     890        sage_free(minimal_cell_reps_of_generators)
     891        sage_free(vertices_to_split)
     892        sage_free(permutation)
     893        sage_free(id_perm)
     894        sage_free(cells_to_refine_by)
     895        sage_free(vertices_determining_current_stack)
     896        PS_dealloc(current_ps)
     897        sage_free(output)
     898        mpz_clear(subgroup_size)
     899        raise MemoryError
     900   
     901    # Copy data of partition to current_ps.
     902    i = 0
     903    j = 0
     904    while partition[i] is not NULL:
     905        k = 0
     906        while partition[i][k] != -1:
     907            current_ps.entries[j+k] = partition[i][k]
     908            current_ps.levels[j+k] = n
     909            k += 1
     910        current_ps.levels[j+k-1] = 0
     911        PS_move_min_to_front(current_ps, j, j+k-1)
     912        j += k
     913        i += 1
     914    current_ps.levels[j-1] = -1
     915    current_ps.depth = 0
     916    current_ps.degree = n
     917   
     918    # default values of "infinity"
     919    for i from 0 <= i < n:
     920        first_indicators[i] = -1
     921    if canonical_label:
     922        for i from 0 <= i < n:
     923            label_indicators[i] = -1
     924   
     925    # set up the identity permutation
     926    for i from 0 <= i < n:
     927        id_perm[i] = i
     928   
     929    # Our first refinement needs to check every cell of the partition,
     930    # so cells_to_refine_by needs to be a list of pointers to each cell.
     931    j = 1
     932    cells_to_refine_by[0] = 0
     933    for i from 0 < i < n:
     934        if current_ps.levels[i-1] == 0:
     935            cells_to_refine_by[j] = i
     936            j += 1
     937    # Ignore the invariant this time, since we are
     938    # creating the root of the search tree.
     939    refine_and_return_invariant(current_ps, S, cells_to_refine_by, j)
     940
     941    # Refine down to a discrete partition
     942    while not PS_is_discrete(current_ps):
     943        if not all_children_are_equivalent(current_ps, S):
     944            current_kids_are_same = current_ps.depth + 1
     945        vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
     946        bitset_unset(vertices_have_been_reduced, current_ps.depth)
     947        i = current_ps.depth
     948        current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     949        first_indicators[i] = current_indicators[i]
     950        if canonical_label:
     951            label_indicators[i] = current_indicators[i]
     952        if base:
     953            output.base[i] = vertices_determining_current_stack[i]
     954
     955    first_meets_current = current_ps.depth
     956    first_kids_are_same = current_ps.depth
     957    first_and_current_indicator_same = current_ps.depth
     958    first_ps = PS_copy(current_ps)
     959    if canonical_label:
     960        label_meets_current = current_ps.depth
     961        label_and_current_indicator_same = current_ps.depth
     962        compared_current_and_label_indicators = 0
     963        label_ps = PS_copy(current_ps)
     964    current_ps.depth -= 1
     965   
     966    # Main loop:
     967    while current_ps.depth != -1:
     968       
     969        # If necessary, update label_meets_current
     970        if canonical_label and label_meets_current > current_ps.depth:
     971            label_meets_current = current_ps.depth
     972
     973        # I. Search for a new vertex to split, and update subgroup information
     974
     975        new_vertex = 0
     976        if current_ps.depth > first_meets_current:
     977            # If we are not at a node of the first stack, reduce size of
     978            # vertices_to_split by using the symmetries we already know.
     979            if not bitset_check(vertices_have_been_reduced, current_ps.depth):
     980                for i from 0 <= i <= index_in_fp_and_mcr:
     981                    j = 0
     982                    while j < current_ps.depth and bitset_check(fixed_points_of_generators[i], vertices_determining_current_stack[j]):
     983                        j += 1
     984                    # If each vertex split so far is fixed by generator i,
     985                    # then remove elements of vertices_to_split which are
     986                    # not minimal in their orbits under generator i.
     987                    if j == current_ps.depth:
     988                        for k from 0 <= k < n:
     989                            if bitset_check(vertices_to_split[current_ps.depth], k) and not bitset_check(minimal_cell_reps_of_generators[i], k):
     990                                bitset_flip(vertices_to_split[current_ps.depth], k)
     991                bitset_flip(vertices_have_been_reduced, current_ps.depth)
     992            # Look for a new point to split.
     993            i = vertices_determining_current_stack[current_ps.depth] + 1
     994            while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
     995                i += 1
     996            if i < n:
     997                # There is a new point.
     998                vertices_determining_current_stack[current_ps.depth] = i
     999                new_vertex = 1
     1000            else:
     1001                # No new point: backtrack.
     1002                current_ps.depth -= 1
     1003        else:
     1004            # If we are at a node of the first stack, the above reduction
     1005            # will not help. Also, we must update information about
     1006            # primary orbits here.
     1007            if current_ps.depth < first_meets_current:
     1008                # If we are done searching under this part of the first
     1009                # stack, then first_meets_current is one higher, and we
     1010                # are looking at a new primary orbit (corresponding to a
     1011                # larger subgroup in the stabilizer chain).
     1012                first_meets_current = current_ps.depth
     1013                for i from 0 <= i < n:
     1014                    if bitset_check(vertices_to_split[current_ps.depth], i):
     1015                        minimal_in_primary_orbit = i
     1016                        break
     1017            while 1:
     1018                i = vertices_determining_current_stack[current_ps.depth]
     1019                # This was the last point to be split here.
     1020                # If it is in the same orbit as minimal_in_primary_orbit,
     1021                # then count it as an element of the primary orbit.
     1022                if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit):
     1023                    subgroup_primary_orbit_size += 1
     1024                # Look for a new point to split.
     1025                i += 1
     1026                while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
     1027                    i += 1
     1028                if i < n:
     1029                    # There is a new point.
     1030                    vertices_determining_current_stack[current_ps.depth] = i
     1031                    if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i:
     1032                        new_vertex = 1
     1033                        break
     1034                else:
     1035                    # No new point: backtrack.
     1036                    # Note that now, we are backtracking up the first stack.
     1037                    vertices_determining_current_stack[current_ps.depth] = -1
     1038                    # If every choice of point to split gave something in the
     1039                    # (same!) primary orbit, then all children of the first
     1040                    # stack at this point are equivalent.
     1041                    j = 0
     1042                    for i from 0 <= i < n:
     1043                        if bitset_check(vertices_to_split[current_ps.depth], i):
     1044                            j += 1
     1045                    if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1:
     1046                        first_kids_are_same = current_ps.depth
     1047                    # Update group size and backtrack.
     1048                    mpz_mul_si(subgroup_size, subgroup_size, subgroup_primary_orbit_size)
     1049                    subgroup_primary_orbit_size = 0
     1050                    current_ps.depth -= 1
     1051                    break
     1052        if not new_vertex:
     1053            continue
     1054
     1055        if current_kids_are_same > current_ps.depth + 1:
     1056            current_kids_are_same = current_ps.depth + 1
     1057        if first_and_current_indicator_same > current_ps.depth:
     1058            first_and_current_indicator_same = current_ps.depth
     1059        if canonical_label and label_and_current_indicator_same > current_ps.depth:
     1060            label_and_current_indicator_same = current_ps.depth
     1061            compared_current_and_label_indicators = 0
     1062       
     1063        # II. Refine down to a discrete partition, or until
     1064        # we leave the part of the tree we are interested in
     1065        discrete = 0
     1066        while 1:
     1067            i = current_ps.depth
     1068            current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     1069            if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]:
     1070                first_and_current_indicator_same += 1
     1071            if canonical_label:
     1072                if compared_current_and_label_indicators == 0:
     1073                    if label_indicators[i] == -1:
     1074                        compared_current_and_label_indicators = -1
     1075                    else:
     1076                        compared_current_and_label_indicators = cmp(current_indicators[i], label_indicators[i])
     1077                if label_and_current_indicator_same == i and compared_current_and_label_indicators == 0:
     1078                    label_and_current_indicator_same += 1
     1079                if compared_current_and_label_indicators > 0:
     1080                    label_indicators[i] = current_indicators[i]
     1081            if first_and_current_indicator_same < current_ps.depth and (not canonical_label or compared_current_and_label_indicators < 0):
     1082                break
     1083            if PS_is_discrete(current_ps):
     1084                discrete = 1
     1085                break
     1086            vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
     1087            bitset_unset(vertices_have_been_reduced, current_ps.depth)
     1088            if not all_children_are_equivalent(current_ps, S):
     1089                current_kids_are_same = current_ps.depth + 1
     1090
     1091        # III. Check for automorphisms and labels
     1092        automorphism = 0
     1093        backtrack = 1 - discrete
     1094        if discrete:
     1095            if current_ps.depth == first_and_current_indicator_same:
     1096                PS_get_perm_from(current_ps, first_ps, permutation)
     1097                if compare_structures(permutation, id_perm, S) == 0:
     1098                    automorphism = 1
     1099        if not automorphism:
     1100            if (not canonical_label) or (compared_current_and_label_indicators < 0):
     1101                backtrack = 1
     1102            else:
     1103                # if we get here, discrete must be true
     1104                update_label = 0
     1105                if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth):
     1106                    update_label = 1
     1107                else:
     1108                    i = compare_structures(current_ps.entries, label_ps.entries, S)
     1109                    if i > 0:
     1110                        update_label = 1
     1111                    elif i < 0:
     1112                        backtrack = 1
     1113                    else:
     1114                        automorphism = 1
     1115                if update_label:
     1116                    PS_copy_from_to(current_ps, label_ps)
     1117                    compared_current_and_label_indicators = 0
     1118                    label_meets_current = current_ps.depth
     1119                    label_and_current_indicator_same = current_ps.depth
     1120                    label_indicators[current_ps.depth] = -1
     1121                    backtrack = 1
     1122        if automorphism:
     1123            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
     1124                index_in_fp_and_mcr += 1
     1125            for i from 0 <= i < n:
     1126                if permutation[i] == i:
     1127                    bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
     1128                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     1129                else:
     1130                    bitset_unset(fixed_points_of_generators[index_in_fp_and_mcr], i)
     1131                    k = i
     1132                    j = permutation[i]
     1133                    while j != i:
     1134                        if j < k: k = j
     1135                        j = permutation[j]
     1136                    if k == i:
     1137                        bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     1138                    else:
     1139                        bitset_unset(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     1140            if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser
     1141                if n*output.num_gens == size_of_generator_array:
     1142                    # must double its size
     1143                    size_of_generator_array *= 2
     1144                    output.generators = <int *> sage_realloc( output.generators, size_of_generator_array )
     1145                    if output.generators is NULL:
     1146                        mem_err = True
     1147                        break # main loop
     1148                j = n*output.num_gens
     1149                for i from 0 <= i < n:
     1150                    output.generators[j + i] = permutation[i]
     1151                output.num_gens += 1
     1152                if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit:
     1153                    current_ps.depth = first_meets_current
     1154                    continue # main loop
     1155            if canonical_label:
     1156                current_ps.depth = label_meets_current
     1157            else:
     1158                current_ps.depth = first_meets_current
     1159            if bitset_check(vertices_have_been_reduced, current_ps.depth):
     1160                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     1161            continue # main loop
     1162        if backtrack:
     1163            i = current_ps.depth
     1164            current_ps.depth = min(max(first_kids_are_same - 1, label_and_current_indicator_same), current_kids_are_same - 1)
     1165            if i == current_kids_are_same:
     1166                continue # main loop
     1167            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
     1168                index_in_fp_and_mcr += 1
     1169            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
     1170            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     1171            j = current_ps.depth
     1172            current_ps.depth = i # just for mcr and fixed functions...
     1173            for i from 0 <= i < n:
     1174                if PS_is_mcr(current_ps, i):
     1175                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     1176                    if PS_is_fixed(current_ps, i):
     1177                        bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
     1178            current_ps.depth = j
     1179            if bitset_check(vertices_have_been_reduced, current_ps.depth):
     1180                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     1181
     1182    # End of main loop.
     1183
     1184    mpz_init_set(output.order, subgroup_size)
     1185    if canonical_label:
     1186        for i from 0 <= i < n:
     1187            output.relabeling[label_ps.entries[i]] = i
     1188   
     1189    # Deallocate:
     1190    for i from 0 <= i < len_of_fp_and_mcr:
     1191        bitset_clear(fixed_points_of_generators[i])
     1192        bitset_clear(minimal_cell_reps_of_generators[i])
     1193    for i from 0 <= i < n:
     1194        bitset_clear(vertices_to_split[i])
     1195    bitset_clear(vertices_have_been_reduced)
     1196    sage_free(current_indicators)
     1197    sage_free(first_indicators)
     1198    if canonical_label:
     1199        PS_dealloc(label_ps)
     1200        sage_free(label_indicators)
     1201    sage_free(fixed_points_of_generators)
     1202    sage_free(minimal_cell_reps_of_generators)
     1203    sage_free(vertices_to_split)
     1204    sage_free(permutation)
     1205    sage_free(id_perm)
     1206    sage_free(cells_to_refine_by)
     1207    sage_free(vertices_determining_current_stack)
     1208    PS_dealloc(current_ps)
     1209    PS_dealloc(first_ps)
     1210    OP_dealloc(orbits_of_subgroup)
     1211    mpz_clear(subgroup_size)
     1212    if not mem_err:
     1213        return output
     1214    else:
     1215        mpz_clear(output.order)
     1216        sage_free(output.generators)
     1217        if base:
     1218            sage_free(output.base)
     1219        if canonical_label:
     1220            sage_free(output.relabeling)
     1221        sage_free(output)
     1222
     1223
     1224
     1225
     1226
     1227
     1228
     1229
     1230
     1231
     1232
     1233
     1234
     1235
     1236
  • setup.py

    diff -r 91af4c3f6b92 -r 074731a26986 setup.py
    a b ext_modules = [ \ 
    675675
    676676    Extension('sage.groups.perm_gps.permgroup_element',
    677677              sources = ['sage/groups/perm_gps/permgroup_element.pyx']), \
     678
     679    Extension('sage.groups.perm_gps.partn_ref.tree_traversal',
     680              sources = ['sage/groups/perm_gps/partn_ref/tree_traversal.pyx']), \
     681
     682    Extension('sage.groups.perm_gps.partn_ref.refinement_graphs',
     683              sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx']), \
    678684
    679685    Extension('sage.structure.sage_object',
    680686              sources = ['sage/structure/sage_object.pyx']), \
    code = setup(name = 'sage', 
    13701376                     'sage.groups.abelian_gps',
    13711377                     'sage.groups.matrix_gps',
    13721378                     'sage.groups.perm_gps',
     1379                     'sage.groups.perm_gps.partn_ref',
    13731380                     
    13741381                     'sage.interfaces',
    13751382