Ticket #3676: trac_3676-graph_isom_refactor.patch

File trac_3676-graph_isom_refactor.patch, 88.0 KB (added by rlm, 13 years ago)

Apply only this patch (please do not delete the others!)

  • 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
    Split up bitset.pxi to allow use in pxd files.
    
    diff -r 5e2426d277eb sage/graphs/graph_generators.py
    a b class DiGraphGenerators(): 
    29912991        import networkx
    29922992        return graph.DiGraph(networkx.gnc_graph(n, seed))
    29932993
     2994    def RandomDirectedGNP(self, n, p):
     2995        r"""
     2996        Returns a random digraph on $n$ nodes.  Each edge is inserted
     2997        independently with probability $p$.
     2998
     2999        REFERENCES:
     3000            [1] P. Erdos and A. Renyi, On Random Graphs, Publ. Math. 6, 290 (1959).
     3001            [2] E. N. Gilbert, Random Graphs, Ann. Math. Stat., 30, 1141 (1959).
     3002       
     3003        PLOTTING:
     3004        When plotting, this graph will use the default spring-layout
     3005        algorithm, unless a position dictionary is specified.
     3006               
     3007        EXAMPLE:
     3008            sage: digraphs.RandomDirectedGNP(10, .2).num_verts()
     3009            10
     3010       
     3011        """
     3012        from random import random
     3013        D = graph.DiGraph(n)
     3014        for i in xrange(n):
     3015            for j in xrange(i):
     3016                if random() < p:
     3017                    D.add_edge(i,j)
     3018            for j in xrange(i+1,n):
     3019                if random() < p:
     3020                    D.add_edge(i,j)
     3021        return D
     3022
    29943023    def RandomDirectedGNR(self, n, p, seed=None):
    29953024        """
    29963025        Returns a random GNR (growing network with redirection) digraph with n
  • new file sage/groups/perm_gps/partn_ref/__init__.py

    diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/__init__.py
    - +  
     1
  • new file sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd

    diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.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'
     11include '../../../misc/bitset_pxd.pxi'
     12
     13from sage.rings.integer cimport Integer
     14
     15cdef struct OrbitPartition:
     16    # Disjoint set data structure for representing the orbits of the generators
     17    # so far found. Also keeps track of the minimum elements of the cells and
     18    # their sizes.
     19    int degree
     20    int *parent
     21    int *rank
     22    int *mcr # minimum cell representatives - only valid at the root of a cell
     23    int *size # also only valid at the root of a cell
     24
     25cdef inline OrbitPartition *OP_new(int)
     26cdef inline int OP_dealloc(OrbitPartition *)
     27cdef inline int OP_find(OrbitPartition *, int)
     28cdef inline int OP_join(OrbitPartition *, int, int)
     29cdef inline int OP_merge_list_perm(OrbitPartition *, int *)
     30
     31cdef struct PartitionStack:
     32    # Representation of a node of the search tree. A sequence of partitions of
     33    # length depth + 1, each of which is finer than the last. Partition k is
     34    # represented as PS.entries in order, broken up immediately after each
     35    # entry of levels which is at most k.
     36    int *entries
     37    int *levels
     38    int depth
     39    int degree
     40
     41cdef inline PartitionStack *PS_new(int, bint)
     42cdef inline PartitionStack *PS_copy(PartitionStack *)
     43cdef inline int PS_copy_from_to(PartitionStack *, PartitionStack *)
     44cdef inline PartitionStack *PS_from_list(object)
     45cdef inline int PS_dealloc(PartitionStack *)
     46cdef inline int PS_print(PartitionStack *)
     47cdef inline int PS_print_partition(PartitionStack *, int)
     48cdef inline bint PS_is_discrete(PartitionStack *)
     49cdef inline int PS_num_cells(PartitionStack *)
     50cdef inline int PS_move_min_to_front(PartitionStack *, int, int)
     51cdef inline bint PS_is_mcr(PartitionStack *, int)
     52cdef inline bint PS_is_fixed(PartitionStack *, int)
     53cdef inline int PS_clear(PartitionStack *)
     54cdef inline int PS_split_point(PartitionStack *, int)
     55cdef inline int PS_first_smallest(PartitionStack *, bitset_t)
     56cdef inline int PS_get_perm_from(PartitionStack *, PartitionStack *, int *)
     57
     58cdef inline int split_point_and_refine(PartitionStack *, int, object,
     59    int (*)(PartitionStack *, object, int *, int), int *)
     60
     61cdef struct aut_gp_and_can_lab_return:
     62    int *generators
     63    int num_gens
     64    int *relabeling
     65    int *base
     66    int base_size
     67    mpz_t order
     68
     69cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab( object, int **, int,
     70    bint (*)(PartitionStack *, object),
     71    int (*)(PartitionStack *, object, int *, int),
     72    int (*)(int *, int *, object), bint, bint, bint)
     73
     74
     75
  • new file sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx

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

    diff -r 5e2426d277eb 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.automorphism_group_canonical_label cimport PartitionStack, PS_is_discrete, PS_move_min_to_front, PS_num_cells, get_aut_gp_and_can_lab, aut_gp_and_can_lab_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 inline int degree(PartitionStack *, CGraph, int, int, bint)
     28cdef inline int sort_by_function(PartitionStack *, int, int *)
     29
     30
     31
  • new file sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

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

    diff -r 5e2426d277eb sage/misc/bitset.pxi
    a b  
    1616
    1717
    1818# This is a .pxi file so that one can inline functions. Doctests in misc_c.
    19 
    20 include "../ext/stdsage.pxi"
    21 
    22 cdef extern from *:
    23     void *memset(void *, int, size_t)
    24     void *memcpy(void *, void *, size_t)
    25     int memcmp(void *, void *, size_t)
    26     size_t strlen(char *)
    27 
    28     # constant literals
    29     int index_shift "(sizeof(unsigned long)==8 ? 6 : 5)"
    30     unsigned long offset_mask "(sizeof(unsigned long)==8 ? 0x3F : 0x1F)"
    31    
    32 cdef struct bitset_s:
    33     long size
    34     long limbs
    35     unsigned long *bits
    36        
    37 ctypedef bitset_s bitset_t[1]
    3819
    3920#############################################################################
    4021# Bitset Initalization
  • new file sage/misc/bitset_pxd.pxi

    diff -r 5e2426d277eb sage/misc/bitset_pxd.pxi
    - +  
     1#*****************************************************************************
     2#     Copyright (C) 2008 Robert Bradshaw <robertwb@math.washington.edu>
     3#
     4#  Distributed under the terms of the GNU General Public License (GPL)
     5#
     6#    This code is distributed in the hope that it will be useful,
     7#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     8#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     9#    General Public License for more details.
     10#
     11#  The full text of the GPL is available at:
     12#
     13#                  http://www.gnu.org/licenses/
     14#*****************************************************************************
     15
     16
     17
     18# This is a .pxi file so that one can inline functions. Doctests in misc_c.
     19
     20include "../ext/stdsage.pxi"
     21
     22cdef extern from *:
     23    void *memset(void *, int, size_t)
     24    void *memcpy(void *, void *, size_t)
     25    int memcmp(void *, void *, size_t)
     26    size_t strlen(char *)
     27
     28    # constant literals
     29    int index_shift "(sizeof(unsigned long)==8 ? 6 : 5)"
     30    unsigned long offset_mask "(sizeof(unsigned long)==8 ? 0x3F : 0x1F)"
     31   
     32cdef struct bitset_s:
     33    long size
     34    long limbs
     35    unsigned long *bits
     36       
     37ctypedef bitset_s bitset_t[1]
     38
  • sage/misc/misc_c.pyx

    diff -r 5e2426d277eb sage/misc/misc_c.pyx
    a b class NonAssociative: 
    277277# Bitset Testing
    278278#############################################################################
    279279
     280include "bitset_pxd.pxi"
    280281include "bitset.pxi"
    281282
    282283def test_bitset(py_a, py_b, long n):
  • setup.py

    diff -r 5e2426d277eb 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.automorphism_group_canonical_label',
     680              sources = ['sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.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