Ticket #4115: trac_4115-double-cosets.patch

File trac_4115-double-cosets.patch, 121.1 KB (added by rlm, 13 years ago)
  • sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd

    # HG changeset patch
    # User Robert L. Miller <rlm@rlmiller.org>
    # Date 1221525292 25200
    # Node ID 2df5cad14396e410bd7abf04853da22ab5f13731
    # Parent  0d158c58d1c23e389772d95e77c543bfa8123446
    Double coset type problems
    
    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd
    a b  
    88
    99include '../../../ext/cdefs.pxi'
    1010include '../../../ext/stdsage.pxi'
    11 include '../../../misc/bitset_pxd.pxi'
     11include 'data_structures_pxd.pxi' # includes bitsets
    1212
    1313from sage.rings.integer cimport Integer
    14 
    15 cdef 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 
    25 cdef inline OrbitPartition *OP_new(int)
    26 cdef inline int OP_dealloc(OrbitPartition *)
    27 cdef inline int OP_find(OrbitPartition *, int)
    28 cdef inline int OP_join(OrbitPartition *, int, int)
    29 cdef inline int OP_merge_list_perm(OrbitPartition *, int *)
    30 
    31 cdef 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 
    41 cdef inline PartitionStack *PS_new(int, bint)
    42 cdef inline PartitionStack *PS_copy(PartitionStack *)
    43 cdef inline int PS_copy_from_to(PartitionStack *, PartitionStack *)
    44 cdef inline PartitionStack *PS_from_list(object)
    45 cdef inline int PS_dealloc(PartitionStack *)
    46 cdef inline int PS_print(PartitionStack *)
    47 cdef inline int PS_print_partition(PartitionStack *, int)
    48 cdef inline bint PS_is_discrete(PartitionStack *)
    49 cdef inline int PS_num_cells(PartitionStack *)
    50 cdef inline int PS_move_min_to_front(PartitionStack *, int, int)
    51 cdef inline bint PS_is_mcr(PartitionStack *, int)
    52 cdef inline bint PS_is_fixed(PartitionStack *, int)
    53 cdef inline int PS_clear(PartitionStack *)
    54 cdef inline int PS_split_point(PartitionStack *, int)
    55 cdef inline int PS_first_smallest(PartitionStack *, bitset_t)
    56 cdef inline int PS_get_perm_from(PartitionStack *, PartitionStack *, int *)
    57 
    58 cdef inline int split_point_and_refine(PartitionStack *, int, object,
    59     int (*)(PartitionStack *, object, int *, int), int *)
    6014
    6115cdef struct aut_gp_and_can_lab_return:
    6216    int *generators
     
    6923cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab( object, int **, int,
    7024    bint (*)(PartitionStack *, object),
    7125    int (*)(PartitionStack *, object, int *, int),
    72     int (*)(int *, int *, object), bint, bint, bint)
     26    int (*)(int *, int *, object, object), bint, bint, bint)
    7327
    7428
    7529
  • sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx
    a b  
    5858   
    5959    Signature:
    6060   
    61     \code{int compare_structures(int *gamma_1, int *gamma_2, object S)}
     61    \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, object S2)}
    6262   
    6363    This function must implement a total ordering on the set of objects of fixed
    64     order. Return -1 if \code{gamma_1(S) < gamma_2(S)}, 0 if
    65     \code{gamma_1(S) == gamma_2(S)}, 1 if \code{gamma_1(S) > gamma_2(S)}.
     64    order. Return -1 if \code{gamma_1(S1) < gamma_2(S2)}, 0 if
     65    \code{gamma_1(S1) == gamma_2(S2)}, 1 if \code{gamma_1(S1) > gamma_2(S2)}.
    6666
    6767C. \code{all_children_are_equivalent}:
    6868   
     
    9393#                         http://www.gnu.org/licenses/
    9494#*****************************************************************************
    9595
    96 include '../../../misc/bitset.pxi'
    97 
    98 cdef inline int min(int a, int b):
    99     if a < b: return a
    100     else: return b
    101 
    102 cdef inline int max(int a, int b):
    103     if a > b: return a
    104     else: return b
     96include 'data_structures_pyx.pxi' # includes bitsets
    10597
    10698cdef inline int cmp(int a, int b):
    10799    if a < b: return -1
    108100    elif a == b: return 0
    109101    else: return 1
    110102
    111 # OrbitPartitions
    112 
    113 cdef inline OrbitPartition *OP_new(int n):
    114     """
    115     Allocate and return a pointer to a new OrbitPartition of degree n. Returns a
    116     null pointer in the case of an allocation failure.
    117     """
    118     cdef int i
    119     cdef OrbitPartition *OP = <OrbitPartition *> \
    120                                 sage_malloc(sizeof(OrbitPartition))
    121     if OP is NULL:
    122         return OP
    123     OP.degree = n
    124     OP.parent = <int *> sage_malloc( n * sizeof(int) )
    125     OP.rank = <int *> sage_malloc( n * sizeof(int) )
    126     OP.mcr = <int *> sage_malloc( n * sizeof(int) )
    127     OP.size = <int *> sage_malloc( n * sizeof(int) )
    128     if OP.parent is NULL or OP.rank is NULL or \
    129        OP.mcr is NULL or OP.size is NULL:
    130         if OP.parent is not NULL: sage_free(OP.parent)
    131         if OP.rank is not NULL: sage_free(OP.rank)
    132         if OP.mcr is not NULL: sage_free(OP.mcr)
    133         if OP.size is not NULL: sage_free(OP.size)
    134         return NULL
    135     for i from 0 <= i < n:
    136         OP.parent[i] = i
    137         OP.rank[i] = 0
    138         OP.mcr[i] = i
    139         OP.size[i] = 1
    140     return OP
    141 
    142 cdef inline int OP_dealloc(OrbitPartition *OP):
    143     sage_free(OP.parent)
    144     sage_free(OP.rank)
    145     sage_free(OP.mcr)
    146     sage_free(OP.size)
    147     sage_free(OP)
    148     return 0
    149 
    150 cdef inline int OP_find(OrbitPartition *OP, int n):
    151     """
    152     Report the representative ("root") of the cell which contains n.
    153     """
    154     if OP.parent[n] == n:
    155         return n
    156     else:
    157         OP.parent[n] = OP_find(OP, OP.parent[n])
    158         return OP.parent[n]
    159 
    160 cdef inline int OP_join(OrbitPartition *OP, int m, int n):
    161     """
    162     Join the cells containing m and n, if they are different.
    163     """
    164     cdef int m_root = OP_find(OP, m)
    165     cdef int n_root = OP_find(OP, n)
    166     if OP.rank[m_root] > OP.rank[n_root]:
    167         OP.parent[n_root] = m_root
    168         OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
    169         OP.size[m_root] += OP.size[n_root]
    170     elif OP.rank[m_root] < OP.rank[n_root]:
    171         OP.parent[m_root] = n_root
    172         OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root])
    173         OP.size[n_root] += OP.size[m_root]
    174     elif m_root != n_root:
    175         OP.parent[n_root] = m_root
    176         OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
    177         OP.size[m_root] += OP.size[n_root]
    178         OP.rank[m_root] += 1
    179 
    180 cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma):
    181     """
    182     Joins the cells of OP which intersect the same orbit of gamma.
    183    
    184     INPUT:
    185         gamma - an integer array representing i -> gamma[i].
    186    
    187     OUTPUT:
    188         1 - something changed
    189         0 - orbits of gamma all contained in cells of OP
    190     """
    191     cdef int i, i_root, gamma_i_root, changed = 0
    192     for i from 0 <= i < OP.degree:
    193         if gamma[i] == i: continue
    194         i_root = OP_find(OP, i)
    195         gamma_i_root = OP_find(OP, gamma[i])
    196         if i_root != gamma_i_root:
    197             changed = 1
    198             OP_join(OP, i_root, gamma_i_root)
    199     return changed
    200 
    201 def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]):
    202     """
    203     Demonstration and testing.
    204    
    205     DOCTEST:
    206         sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
    207         ... import OP_represent
    208         sage: OP_represent()
    209         Allocating OrbitPartition...
    210         Allocation passed.
    211         Checking that each element reports itself as its root.
    212         Each element reports itself as its root.
    213         Merging:
    214         Merged 0 and 1.
    215         Merged 2 and 3.
    216         Merged 3 and 4.
    217         Done merging.
    218         Finding:
    219         0 -> 0, root: size=2, mcr=0, rank=1
    220         1 -> 0
    221         2 -> 2, root: size=3, mcr=2, rank=1
    222         3 -> 2
    223         4 -> 2
    224         5 -> 5, root: size=1, mcr=5, rank=0
    225         6 -> 6, root: size=1, mcr=6, rank=0
    226         7 -> 7, root: size=1, mcr=7, rank=0
    227         8 -> 8, root: size=1, mcr=8, rank=0
    228         Allocating array to test merge_perm.
    229         Allocation passed.
    230         Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8]
    231         Done merging.
    232         Finding:
    233         0 -> 0, root: size=5, mcr=0, rank=2
    234         1 -> 0
    235         2 -> 0
    236         3 -> 0
    237         4 -> 0
    238         5 -> 5, root: size=3, mcr=5, rank=1
    239         6 -> 5
    240         7 -> 5
    241         8 -> 8, root: size=1, mcr=8, rank=0
    242         Deallocating OrbitPartition.
    243         Done.
    244 
    245     """
    246     cdef int i
    247     print "Allocating OrbitPartition..."
    248     cdef OrbitPartition *OP = OP_new(n)
    249     if OP is NULL:
    250         print "Allocation failed!"
    251         return
    252     print "Allocation passed."
    253     print "Checking that each element reports itself as its root."
    254     good = True
    255     for i from 0 <= i < n:
    256         if not OP_find(OP, i) == i:
    257             print "Failed at i = %d!"%i
    258             good = False
    259     if good: print "Each element reports itself as its root."
    260     print "Merging:"
    261     for i,j in merges:
    262         OP_join(OP, i, j)
    263         print "Merged %d and %d."%(i,j)
    264     print "Done merging."
    265     print "Finding:"
    266     for i from 0 <= i < n:
    267         j = OP_find(OP, i)
    268         s = "%d -> %d"%(i, j)
    269         if i == j:
    270             s += ", root: size=%d, mcr=%d, rank=%d"%\
    271                    (OP.size[i], OP.mcr[i], OP.rank[i])
    272         print s
    273     print "Allocating array to test merge_perm."
    274     cdef int *gamma = <int *> sage_malloc( n * sizeof(int) )
    275     if gamma is NULL:
    276         print "Allocation failed!"
    277         OP_dealloc(OP)
    278         return
    279     print "Allocation passed."
    280     for i from 0 <= i < n:
    281         gamma[i] = perm[i]
    282     print "Merging permutation: %s"%perm
    283     OP_merge_list_perm(OP, gamma)
    284     print "Done merging."
    285     print "Finding:"
    286     for i from 0 <= i < n:
    287         j = OP_find(OP, i)
    288         s = "%d -> %d"%(i, j)
    289         if i == j:
    290             s += ", root: size=%d, mcr=%d, rank=%d"%\
    291                    (OP.size[i], OP.mcr[i], OP.rank[i])
    292         print s   
    293     print "Deallocating OrbitPartition."
    294     sage_free(gamma)
    295     OP_dealloc(OP)
    296     print "Done."
    297 
    298 # PartitionStacks
    299 
    300 cdef inline PartitionStack *PS_new(int n, bint unit_partition):
    301     """
    302     Allocate and return a pointer to a new PartitionStack of degree n. Returns a
    303     null pointer in the case of an allocation failure.
    304     """
    305     cdef int i
    306     cdef PartitionStack *PS = <PartitionStack *> \
    307                                 sage_malloc(sizeof(PartitionStack))
    308     if PS is NULL:
    309         return PS
    310     PS.entries = <int *> sage_malloc( n * sizeof(int) )
    311     PS.levels = <int *> sage_malloc( n * sizeof(int) )
    312     if PS.entries is NULL or PS.levels is NULL:
    313         if PS.entries is not NULL: sage_free(PS.entries)
    314         if PS.levels is not NULL: sage_free(PS.levels)
    315         return NULL
    316     PS.depth = 0
    317     PS.degree = n
    318     if unit_partition:
    319         for i from 0 <= i < n-1:
    320             PS.entries[i] = i
    321             PS.levels[i] = n
    322         PS.entries[n-1] = n-1
    323         PS.levels[n-1] = -1
    324     return PS
    325 
    326 cdef inline PartitionStack *PS_copy(PartitionStack *PS):
    327     """
    328     Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
    329     pointer in the case of an allocation failure.
    330     """
    331     cdef int i
    332     cdef PartitionStack *PS2 = <PartitionStack *> \
    333                                 sage_malloc(sizeof(PartitionStack))
    334     if PS2 is NULL:
    335         return PS2
    336     PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
    337     PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
    338     if PS2.entries is NULL or PS2.levels is NULL:
    339         if PS2.entries is not NULL: sage_free(PS2.entries)
    340         if PS2.levels is not NULL: sage_free(PS2.levels)
    341         return NULL
    342     PS_copy_from_to(PS, PS2)
    343     return PS2
    344 
    345 cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2):
    346     """
    347     Copy all data from PS to PS2.
    348     """
    349     PS2.depth = PS.depth
    350     PS2.degree = PS.degree
    351     for i from 0 <= i < PS.degree:
    352         PS2.entries[i] = PS.entries[i]
    353         PS2.levels[i] = PS.levels[i]   
    354 
    355 cdef inline PartitionStack *PS_from_list(object L):
    356     """
    357     Allocate and return a pointer to a PartitionStack representing L. Returns a
    358     null pointer in the case of an allocation failure.
    359     """
    360     cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0
    361     for cell from 0 <= cell < num_cells:
    362         n += len(L[cell])
    363     cdef PartitionStack *PS = PS_new(n, 0)
    364     cell = 0
    365     if PS is NULL:
    366         return PS
    367     while 1:
    368         cur_len = len(L[cell])
    369         for i from 0 <= i < cur_len:
    370             PS.entries[cur_start + i] = L[cell][i]
    371             PS.levels[cur_start + i] = n
    372         PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
    373         cur_start += cur_len
    374         cell += 1
    375         if cell == num_cells:
    376             PS.levels[cur_start-1] = -1
    377             break
    378         PS.levels[cur_start-1] = 0
    379     return PS
    380 
    381 cdef inline int PS_dealloc(PartitionStack *PS):
    382     sage_free(PS.entries)
    383     sage_free(PS.levels)
    384     sage_free(PS)
    385     return 0
    386 
    387 cdef inline int PS_print(PartitionStack *PS):
    388     """
    389     Print a visual representation of PS.
    390     """
    391     cdef int i
    392     for i from 0 <= i < PS.depth + 1:
    393         PS_print_partition(PS, i)
    394     print
    395 
    396 cdef inline int PS_print_partition(PartitionStack *PS, int k):
    397     """
    398     Print the partition at depth k.
    399     """
    400     s = '('
    401     for i from 0 <= i < PS.degree:
    402         s += str(PS.entries[i])
    403         if PS.levels[i] <= k:
    404             s += '|'
    405         else:
    406             s += ','
    407     s = s[:-1] + ')'
    408     print s
    409 
    410 cdef inline bint PS_is_discrete(PartitionStack *PS):
    411     """
    412     Returns whether the deepest partition consists only of singleton cells.
    413     """
    414     cdef int i
    415     for i from 0 <= i < PS.degree:
    416         if PS.levels[i] > PS.depth:
    417             return 0
    418     return 1
    419 
    420 cdef inline int PS_num_cells(PartitionStack *PS):
    421     """
    422     Returns the number of cells.
    423     """
    424     cdef int i, ncells = 0
    425     for i from 0 <= i < PS.degree:
    426         if PS.levels[i] <= PS.depth:
    427             ncells += 1
    428     return ncells
    429 
    430 cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
    431     """
    432     Makes sure that the first element of the segment of entries i with
    433     start <= i <= end is minimal.
    434     """
    435     cdef int i, min_loc = start, minimum = PS.entries[start]
    436     for i from start < i <= end:
    437         if PS.entries[i] < minimum:
    438             min_loc = i
    439             minimum = PS.entries[i]
    440     if min_loc != start:
    441         PS.entries[min_loc] = PS.entries[start]
    442         PS.entries[start] = minimum
    443 
    444 cdef inline bint PS_is_mcr(PartitionStack *PS, int m):
    445     """
    446     Returns whether PS.elements[m] (not m!) is the smallest element of its cell.
    447     """
    448     return m == 0 or PS.levels[m-1] <= PS.depth
    449 
    450 cdef inline bint PS_is_fixed(PartitionStack *PS, int m):
    451     """
    452     Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming
    453     PS_is_mcr(PS, m) is already true.
    454     """
    455     return PS.levels[m] <= PS.depth
    456 
    457 cdef inline int PS_clear(PartitionStack *PS):
    458     """
    459     Sets the current partition to the first shallower one, i.e. forgets about
    460     boundaries between cells that are new to the current level.
    461     """
    462     cdef int i, cur_start = 0
    463     for i from 0 <= i < PS.degree:
    464         if PS.levels[i] >= PS.depth:
    465             PS.levels[i] += 1
    466         if PS.levels[i] < PS.depth:
    467             PS_move_min_to_front(PS, cur_start, i)
    468             cur_start = i+1
    469 
    470 cdef inline int PS_split_point(PartitionStack *PS, int v):
    471     """
    472     Detaches the point v from the cell it is in, putting the singleton cell
    473     of just v in front. Returns the position where v is now located.
    474     """
    475     cdef int i = 0, index_of_v
    476     while PS.entries[i] != v:
    477         i += 1
    478     index_of_v = i
    479     while PS.levels[i] > PS.depth:
    480         i += 1
    481     if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \
    482        and PS.levels[index_of_v] > PS.depth:
    483         # if v is first (make sure v is not already alone)
    484         PS_move_min_to_front(PS, index_of_v+1, i)
    485         PS.levels[index_of_v] = PS.depth
    486         return index_of_v
    487     else:
    488         # If v is not at front, i.e. v is not minimal in its cell,
    489         # then move_min_to_front is not necessary since v will swap
    490         # with the first before making its own cell, leaving it at
    491         # the front of the other.
    492         i = index_of_v
    493         while i != 0 and PS.levels[i-1] > PS.depth:
    494             i -= 1
    495         PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v
    496         PS.entries[i+1] = PS.entries[i] # move the first (min) to second
    497         PS.entries[i] = v # place v first
    498         PS.levels[i] = PS.depth
    499         return i
    500 
    501 cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
    502     """
    503     Find the first occurrence of the smallest cell of size greater than one,
    504     store its entries to b, and return its minimum element.
    505     """
    506     cdef int i = 0, j = 0, location = 0, n = PS.degree
    507     bitset_zero(b)
    508     while 1:
    509         if PS.levels[i] <= PS.depth:
    510             if i != j and n > i - j + 1:
    511                 n = i - j + 1
    512                 location = j
    513             j = i + 1
    514         if PS.levels[i] == -1: break
    515         i += 1
    516     # location now points to the beginning of the first, smallest,
    517     # nontrivial cell
    518     i = location
    519     while 1:
    520         bitset_flip(b, PS.entries[i])
    521         if PS.levels[i] <= PS.depth: break
    522         i += 1
    523     return PS.entries[location]
    524 
    525 cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
    526     """
    527     Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
    528     denotes the entry of the ith cell of the discrete partition PS.
    529     """
    530     cdef int i = 0
    531     for i from 0 <= i < PS1.degree:
    532         gamma[PS2.entries[i]] = PS1.entries[i]
    533 
    534 def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]):
    535     """
    536     Demonstration and testing.
    537    
    538     DOCTEST:
    539         sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
    540         ... import PS_represent
    541         sage: PS_represent()
    542         Allocating PartitionStack...
    543         Allocation passed:
    544         (0,1,2,3,4,5,6,7,8,9)
    545         Checking that entries are in order and correct level.
    546         Everything seems in order, deallocating.
    547         Deallocated.
    548         Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]].
    549         PartitionStack's data:
    550         entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2]
    551         levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1]
    552         depth = 0, degree = 10
    553         (6|3,4,8,7|1,9,5|0,2)
    554         Checking PS_is_discrete:
    555         False
    556         Checking PS_num_cells:
    557         4
    558         Checking PS_is_mcr, min cell reps are:
    559         [6, 3, 1, 0]
    560         Checking PS_is_fixed, fixed elements are:
    561         [6]
    562         Copying PartitionStack:
    563         (6|3,4,8,7|1,9,5|0,2)
    564         Checking for consistency.
    565         Everything is consistent.
    566         Clearing copy:
    567         (0,3,4,8,7,1,9,5,6,2)
    568         Splitting point 6 from original:
    569         0
    570         (6|3,4,8,7|1,9,5|0,2)
    571         Splitting point 1 from original:
    572         5
    573         (6|3,4,8,7|1|5,9|0,2)
    574         Splitting point 8 from original:
    575         1
    576         (6|8|3,4,7|1|5,9|0,2)
    577         Splitting point 2 from original:
    578         8
    579         (6|8|3,4,7|1|5,9|2|0)
    580         Getting permutation from PS2->PS:
    581         [6, 1, 0, 8, 3, 9, 2, 7, 4, 5]
    582         Finding first smallest:
    583         Minimal element is 5, bitset is:
    584         0000010001
    585         Deallocating PartitionStacks.
    586         Done.
    587 
    588     """
    589     cdef int i, n = sum([len(cell) for cell in partition]), *gamma
    590     cdef bitset_t b
    591     print "Allocating PartitionStack..."
    592     cdef PartitionStack *PS = PS_new(n, 1), *PS2
    593     if PS is NULL:
    594         print "Allocation failed!"
    595         return
    596     print "Allocation passed:"
    597     PS_print(PS)
    598     print "Checking that entries are in order and correct level."
    599     good = True
    600     for i from 0 <= i < n-1:
    601         if not (PS.entries[i] == i and PS.levels[i] == n):
    602             print "Failed at i = %d!"%i
    603             print PS.entries[i], PS.levels[i], i, n
    604             good = False
    605     if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1):
    606         print "Failed at i = %d!"%(n-1)
    607         good = False
    608     if not PS.degree == n or not PS.depth == 0:
    609         print "Incorrect degree or depth!"
    610         good = False
    611     if good: print "Everything seems in order, deallocating."
    612     PS_dealloc(PS)
    613     print "Deallocated."
    614     print "Creating PartitionStack from partition %s."%partition
    615     PS = PS_from_list(partition)
    616     print "PartitionStack's data:"
    617     print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n]
    618     print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n]
    619     print "depth = %d, degree = %d"%(PS.depth,PS.degree)
    620     PS_print(PS)
    621     print "Checking PS_is_discrete:"
    622     print "True" if PS_is_discrete(PS) else "False"
    623     print "Checking PS_num_cells:"
    624     print PS_num_cells(PS)
    625     print "Checking PS_is_mcr, min cell reps are:"
    626     L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)]
    627     print L
    628     print "Checking PS_is_fixed, fixed elements are:"
    629     print [PS.entries[l] for l in L if PS_is_fixed(PS, l)]
    630     print "Copying PartitionStack:"
    631     PS2 = PS_copy(PS)
    632     PS_print(PS2)
    633     print "Checking for consistency."
    634     good = True
    635     for i from 0 <= i < n:
    636         if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]:
    637             print "Failed at i = %d!"%i
    638             good = False
    639     if PS.degree != PS2.degree or PS.depth != PS2.depth:
    640         print "Failure with degree or depth!"
    641         good = False
    642     if good:
    643         print "Everything is consistent."
    644     print "Clearing copy:"
    645     PS_clear(PS2)
    646     PS_print(PS2)
    647     for s in splits:
    648         print "Splitting point %d from original:"%s
    649         print PS_split_point(PS, s)
    650         PS_print(PS)
    651     print "Getting permutation from PS2->PS:"
    652     gamma = <int *> sage_malloc(n * sizeof(int))
    653     PS_get_perm_from(PS, PS2, gamma)
    654     print [gamma[i] for i from 0 <= i < n]
    655     sage_free(gamma)
    656     print "Finding first smallest:"
    657     bitset_init(b, n)
    658     i = PS_first_smallest(PS, b)
    659     print "Minimal element is %d, bitset is:"%i
    660     print bitset_string(b)
    661     bitset_clear(b)
    662     print "Deallocating PartitionStacks."
    663     PS_dealloc(PS)
    664     PS_dealloc(PS2)
    665     print "Done."
    666 
    667103# Functions
    668 
    669 cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
    670     int (*refine_and_return_invariant)\
    671          (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
    672     int *cells_to_refine_by):
    673     """
    674     Make the partition stack one longer by copying the last partition in the
    675     stack, split off a given point, and refine. Return the invariant given by
    676     the refinement function.
    677    
    678     INPUT:
    679     PS -- the partition stack to refine
    680     v -- the point to split
    681     S -- the structure
    682     refine_and_return_invariant -- the refinement function provided
    683     cells_to_refine_by -- an array, contents ignored
    684 
    685     """
    686     PS.depth += 1
    687     PS_clear(PS)
    688     cells_to_refine_by[0] = PS_split_point(PS, v)
    689     return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
    690104
    691105cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, object S):
    692106    return 0
     
    694108cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
    695109    return 0
    696110
    697 cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S):
     111cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S1, object S2):
    698112    return 0
    699113
    700114def test_get_aut_gp_and_can_lab_trivially(int n=6, partition=[[0,1,2],[3,4],[5]],
     
    747161    bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
    748162    int (*refine_and_return_invariant)\
    749163         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
    750     int (*compare_structures)(int *gamma_1, int *gamma_2, object S),
     164    int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2),
    751165    bint canonical_label, bint base, bint order):
    752166    """
    753167    Traverse the search space for subgroup/canonical label calculation.
     
    779193        int -- returns an invariant under application of arbitrary permutations
    780194    compare_structures -- pointer to a function
    781195        INPUT:
    782         gamma_1, gamma_2 -- (list) permutations of the points of S
    783         S -- pointer to the structure
     196        gamma_1, gamma_2 -- (list) permutations of the points of S1 and S2
     197        S1, S2 -- pointers to the structures
    784198        OUTPUT:
    785         int -- 0 if gamma_1(S) = gamma_2(S), otherwise -1 or 1 (see docs for cmp),
     199        int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp),
    786200            such that the set of all structures is well-ordered
    787201    OUTPUT:
    788202    pointer to a aut_gp_and_can_lab_return struct
     
    1020434    # Ignore the invariant this time, since we are
    1021435    # creating the root of the search tree.
    1022436    refine_and_return_invariant(current_ps, S, cells_to_refine_by, j)
    1023 
     437    PS_move_all_mins_to_front(current_ps)
     438   
    1024439    # Refine down to a discrete partition
    1025440    while not PS_is_discrete(current_ps):
    1026441        if not all_children_are_equivalent(current_ps, S):
     
    1029444        bitset_unset(vertices_have_been_reduced, current_ps.depth)
    1030445        i = current_ps.depth
    1031446        current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     447        PS_move_all_mins_to_front(current_ps)
    1032448        first_indicators[i] = current_indicators[i]
    1033449        if canonical_label:
    1034450            label_indicators[i] = current_indicators[i]
     
    1074490                bitset_flip(vertices_have_been_reduced, current_ps.depth)
    1075491            # Look for a new point to split.
    1076492            i = vertices_determining_current_stack[current_ps.depth] + 1
    1077             while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
    1078                 i += 1
    1079             if i < n:
     493            i = bitset_next(vertices_to_split[current_ps.depth], i)
     494            if i != -1:
    1080495                # There is a new point.
    1081496                vertices_determining_current_stack[current_ps.depth] = i
    1082497                new_vertex = 1
     
    1106521                    subgroup_primary_orbit_size += 1
    1107522                # Look for a new point to split.
    1108523                i += 1
    1109                 while i < n and not bitset_check(vertices_to_split[current_ps.depth], i):
    1110                     i += 1
    1111                 if i < n:
     524                i = bitset_next(vertices_to_split[current_ps.depth], i)
     525                if i != -1:
    1112526                    # There is a new point.
    1113527                    vertices_determining_current_stack[current_ps.depth] = i
    1114528                    if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i:
     
    1139553            current_kids_are_same = current_ps.depth + 1
    1140554        if first_and_current_indicator_same > current_ps.depth:
    1141555            first_and_current_indicator_same = current_ps.depth
    1142         if canonical_label and label_and_current_indicator_same > current_ps.depth:
     556        if canonical_label and label_and_current_indicator_same >= current_ps.depth:
    1143557            label_and_current_indicator_same = current_ps.depth
    1144558            compared_current_and_label_indicators = 0
    1145559       
     
    1149563        while 1:
    1150564            i = current_ps.depth
    1151565            current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by)
     566            PS_move_all_mins_to_front(current_ps)
    1152567            if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]:
    1153568                first_and_current_indicator_same += 1
    1154569            if canonical_label:
     
    1177592        if discrete:
    1178593            if current_ps.depth == first_and_current_indicator_same:
    1179594                PS_get_perm_from(current_ps, first_ps, permutation)
    1180                 if compare_structures(permutation, id_perm, S) == 0:
     595                if compare_structures(permutation, id_perm, S, S) == 0:
    1181596                    automorphism = 1
    1182597        if not automorphism:
    1183598            if (not canonical_label) or (compared_current_and_label_indicators < 0):
     
    1188603                if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth):
    1189604                    update_label = 1
    1190605                else:
    1191                     i = compare_structures(current_ps.entries, label_ps.entries, S)
     606                    i = compare_structures(current_ps.entries, label_ps.entries, S, S)
    1192607                    if i > 0:
    1193608                        update_label = 1
    1194609                    elif i < 0:
  • new file sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi
    - +  
     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
     13cdef struct OrbitPartition:
     14    # Disjoint set data structure for representing the orbits of the generators
     15    # so far found. Also keeps track of the minimum elements of the cells and
     16    # their sizes.
     17    int degree
     18    int *parent
     19    int *rank
     20    int *mcr # minimum cell representatives - only valid at the root of a cell
     21    int *size # also only valid at the root of a cell
     22
     23cdef struct PartitionStack:
     24    # Representation of a node of the search tree. A sequence of partitions of
     25    # length depth + 1, each of which is finer than the last. Partition k is
     26    # represented as PS.entries in order, broken up immediately after each
     27    # entry of levels which is at most k.
     28    int *entries
     29    int *levels
     30    int depth
     31    int degree
     32
     33
  • new file sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi
    - +  
     1r"""
     2Data structures
     3
     4This module implements TODO...
     5
     6REFERENCE:
     7
     8    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
     9        Vol. 30 (1981), pp. 45-87.
     10
     11"""
     12
     13#*****************************************************************************
     14#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     15#
     16# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     17#                         http://www.gnu.org/licenses/
     18#*****************************************************************************
     19
     20include '../../../misc/bitset.pxi'
     21
     22cdef inline int min(int a, int b):
     23    if a < b: return a
     24    else: return b
     25
     26cdef inline int max(int a, int b):
     27    if a > b: return a
     28    else: return b
     29
     30# OrbitPartitions
     31
     32cdef inline OrbitPartition *OP_new(int n):
     33    """
     34    Allocate and return a pointer to a new OrbitPartition of degree n. Returns a
     35    null pointer in the case of an allocation failure.
     36    """
     37    cdef int i
     38    cdef OrbitPartition *OP = <OrbitPartition *> \
     39                                sage_malloc(sizeof(OrbitPartition))
     40    if OP is NULL:
     41        return OP
     42    OP.degree = n
     43    OP.parent = <int *> sage_malloc( n * sizeof(int) )
     44    OP.rank = <int *> sage_malloc( n * sizeof(int) )
     45    OP.mcr = <int *> sage_malloc( n * sizeof(int) )
     46    OP.size = <int *> sage_malloc( n * sizeof(int) )
     47    if OP.parent is NULL or OP.rank is NULL or \
     48       OP.mcr is NULL or OP.size is NULL:
     49        if OP.parent is not NULL: sage_free(OP.parent)
     50        if OP.rank is not NULL: sage_free(OP.rank)
     51        if OP.mcr is not NULL: sage_free(OP.mcr)
     52        if OP.size is not NULL: sage_free(OP.size)
     53        return NULL
     54    for i from 0 <= i < n:
     55        OP.parent[i] = i
     56        OP.rank[i] = 0
     57        OP.mcr[i] = i
     58        OP.size[i] = 1
     59    return OP
     60
     61cdef inline int OP_dealloc(OrbitPartition *OP):
     62    sage_free(OP.parent)
     63    sage_free(OP.rank)
     64    sage_free(OP.mcr)
     65    sage_free(OP.size)
     66    sage_free(OP)
     67    return 0
     68
     69cdef inline int OP_find(OrbitPartition *OP, int n):
     70    """
     71    Report the representative ("root") of the cell which contains n.
     72    """
     73    if OP.parent[n] == n:
     74        return n
     75    else:
     76        OP.parent[n] = OP_find(OP, OP.parent[n])
     77        return OP.parent[n]
     78
     79cdef inline int OP_join(OrbitPartition *OP, int m, int n):
     80    """
     81    Join the cells containing m and n, if they are different.
     82    """
     83    cdef int m_root = OP_find(OP, m)
     84    cdef int n_root = OP_find(OP, n)
     85    if OP.rank[m_root] > OP.rank[n_root]:
     86        OP.parent[n_root] = m_root
     87        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
     88        OP.size[m_root] += OP.size[n_root]
     89    elif OP.rank[m_root] < OP.rank[n_root]:
     90        OP.parent[m_root] = n_root
     91        OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root])
     92        OP.size[n_root] += OP.size[m_root]
     93    elif m_root != n_root:
     94        OP.parent[n_root] = m_root
     95        OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root])
     96        OP.size[m_root] += OP.size[n_root]
     97        OP.rank[m_root] += 1
     98
     99cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma):
     100    """
     101    Joins the cells of OP which intersect the same orbit of gamma.
     102   
     103    INPUT:
     104        gamma - an integer array representing i -> gamma[i].
     105   
     106    OUTPUT:
     107        1 - something changed
     108        0 - orbits of gamma all contained in cells of OP
     109    """
     110    cdef int i, i_root, gamma_i_root, changed = 0
     111    for i from 0 <= i < OP.degree:
     112        if gamma[i] == i: continue
     113        i_root = OP_find(OP, i)
     114        gamma_i_root = OP_find(OP, gamma[i])
     115        if i_root != gamma_i_root:
     116            changed = 1
     117            OP_join(OP, i_root, gamma_i_root)
     118    return changed
     119
     120def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]):
     121    """
     122    Demonstration and testing.
     123   
     124    DOCTEST:
     125        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
     126        ... import OP_represent
     127        sage: OP_represent()
     128        Allocating OrbitPartition...
     129        Allocation passed.
     130        Checking that each element reports itself as its root.
     131        Each element reports itself as its root.
     132        Merging:
     133        Merged 0 and 1.
     134        Merged 2 and 3.
     135        Merged 3 and 4.
     136        Done merging.
     137        Finding:
     138        0 -> 0, root: size=2, mcr=0, rank=1
     139        1 -> 0
     140        2 -> 2, root: size=3, mcr=2, rank=1
     141        3 -> 2
     142        4 -> 2
     143        5 -> 5, root: size=1, mcr=5, rank=0
     144        6 -> 6, root: size=1, mcr=6, rank=0
     145        7 -> 7, root: size=1, mcr=7, rank=0
     146        8 -> 8, root: size=1, mcr=8, rank=0
     147        Allocating array to test merge_perm.
     148        Allocation passed.
     149        Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8]
     150        Done merging.
     151        Finding:
     152        0 -> 0, root: size=5, mcr=0, rank=2
     153        1 -> 0
     154        2 -> 0
     155        3 -> 0
     156        4 -> 0
     157        5 -> 5, root: size=3, mcr=5, rank=1
     158        6 -> 5
     159        7 -> 5
     160        8 -> 8, root: size=1, mcr=8, rank=0
     161        Deallocating OrbitPartition.
     162        Done.
     163
     164    """
     165    cdef int i
     166    print "Allocating OrbitPartition..."
     167    cdef OrbitPartition *OP = OP_new(n)
     168    if OP is NULL:
     169        print "Allocation failed!"
     170        return
     171    print "Allocation passed."
     172    print "Checking that each element reports itself as its root."
     173    good = True
     174    for i from 0 <= i < n:
     175        if not OP_find(OP, i) == i:
     176            print "Failed at i = %d!"%i
     177            good = False
     178    if good: print "Each element reports itself as its root."
     179    print "Merging:"
     180    for i,j in merges:
     181        OP_join(OP, i, j)
     182        print "Merged %d and %d."%(i,j)
     183    print "Done merging."
     184    print "Finding:"
     185    for i from 0 <= i < n:
     186        j = OP_find(OP, i)
     187        s = "%d -> %d"%(i, j)
     188        if i == j:
     189            s += ", root: size=%d, mcr=%d, rank=%d"%\
     190                   (OP.size[i], OP.mcr[i], OP.rank[i])
     191        print s
     192    print "Allocating array to test merge_perm."
     193    cdef int *gamma = <int *> sage_malloc( n * sizeof(int) )
     194    if gamma is NULL:
     195        print "Allocation failed!"
     196        OP_dealloc(OP)
     197        return
     198    print "Allocation passed."
     199    for i from 0 <= i < n:
     200        gamma[i] = perm[i]
     201    print "Merging permutation: %s"%perm
     202    OP_merge_list_perm(OP, gamma)
     203    print "Done merging."
     204    print "Finding:"
     205    for i from 0 <= i < n:
     206        j = OP_find(OP, i)
     207        s = "%d -> %d"%(i, j)
     208        if i == j:
     209            s += ", root: size=%d, mcr=%d, rank=%d"%\
     210                   (OP.size[i], OP.mcr[i], OP.rank[i])
     211        print s   
     212    print "Deallocating OrbitPartition."
     213    sage_free(gamma)
     214    OP_dealloc(OP)
     215    print "Done."
     216
     217# PartitionStacks
     218
     219cdef inline PartitionStack *PS_new(int n, bint unit_partition):
     220    """
     221    Allocate and return a pointer to a new PartitionStack of degree n. Returns a
     222    null pointer in the case of an allocation failure.
     223    """
     224    cdef int i
     225    cdef PartitionStack *PS = <PartitionStack *> \
     226                                sage_malloc(sizeof(PartitionStack))
     227    if PS is NULL:
     228        return PS
     229    PS.entries = <int *> sage_malloc( n * sizeof(int) )
     230    PS.levels = <int *> sage_malloc( n * sizeof(int) )
     231    if PS.entries is NULL or PS.levels is NULL:
     232        if PS.entries is not NULL: sage_free(PS.entries)
     233        if PS.levels is not NULL: sage_free(PS.levels)
     234        return NULL
     235    PS.depth = 0
     236    PS.degree = n
     237    if unit_partition:
     238        for i from 0 <= i < n-1:
     239            PS.entries[i] = i
     240            PS.levels[i] = n
     241        PS.entries[n-1] = n-1
     242        PS.levels[n-1] = -1
     243    return PS
     244
     245cdef inline PartitionStack *PS_copy(PartitionStack *PS):
     246    """
     247    Allocate and return a pointer to a copy of PartitionStack PS. Returns a null
     248    pointer in the case of an allocation failure.
     249    """
     250    cdef int i
     251    cdef PartitionStack *PS2 = <PartitionStack *> \
     252                                sage_malloc(sizeof(PartitionStack))
     253    if PS2 is NULL:
     254        return PS2
     255    PS2.entries = <int *> sage_malloc( PS.degree * sizeof(int) )
     256    PS2.levels = <int *> sage_malloc( PS.degree * sizeof(int) )
     257    if PS2.entries is NULL or PS2.levels is NULL:
     258        if PS2.entries is not NULL: sage_free(PS2.entries)
     259        if PS2.levels is not NULL: sage_free(PS2.levels)
     260        return NULL
     261    PS_copy_from_to(PS, PS2)
     262    return PS2
     263
     264cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2):
     265    """
     266    Copy all data from PS to PS2.
     267    """
     268    PS2.depth = PS.depth
     269    PS2.degree = PS.degree
     270    for i from 0 <= i < PS.degree:
     271        PS2.entries[i] = PS.entries[i]
     272        PS2.levels[i] = PS.levels[i]   
     273
     274cdef inline PartitionStack *PS_from_list(object L):
     275    """
     276    Allocate and return a pointer to a PartitionStack representing L. Returns a
     277    null pointer in the case of an allocation failure.
     278    """
     279    cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0
     280    for cell from 0 <= cell < num_cells:
     281        n += len(L[cell])
     282    cdef PartitionStack *PS = PS_new(n, 0)
     283    cell = 0
     284    if PS is NULL:
     285        return PS
     286    while 1:
     287        cur_len = len(L[cell])
     288        for i from 0 <= i < cur_len:
     289            PS.entries[cur_start + i] = L[cell][i]
     290            PS.levels[cur_start + i] = n
     291        PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1)
     292        cur_start += cur_len
     293        cell += 1
     294        if cell == num_cells:
     295            PS.levels[cur_start-1] = -1
     296            break
     297        PS.levels[cur_start-1] = 0
     298    return PS
     299
     300cdef inline int PS_dealloc(PartitionStack *PS):
     301    sage_free(PS.entries)
     302    sage_free(PS.levels)
     303    sage_free(PS)
     304    return 0
     305
     306cdef inline int PS_print(PartitionStack *PS):
     307    """
     308    Print a visual representation of PS.
     309    """
     310    cdef int i
     311    for i from 0 <= i < PS.depth + 1:
     312        PS_print_partition(PS, i)
     313    print
     314
     315cdef inline int PS_print_partition(PartitionStack *PS, int k):
     316    """
     317    Print the partition at depth k.
     318    """
     319    s = '('
     320    for i from 0 <= i < PS.degree:
     321        s += str(PS.entries[i])
     322        if PS.levels[i] <= k:
     323            s += '|'
     324        else:
     325            s += ','
     326    s = s[:-1] + ')'
     327    print s
     328
     329cdef inline bint PS_is_discrete(PartitionStack *PS):
     330    """
     331    Returns whether the deepest partition consists only of singleton cells.
     332    """
     333    cdef int i
     334    for i from 0 <= i < PS.degree:
     335        if PS.levels[i] > PS.depth:
     336            return 0
     337    return 1
     338
     339cdef inline int PS_num_cells(PartitionStack *PS):
     340    """
     341    Returns the number of cells.
     342    """
     343    cdef int i, ncells = 0
     344    for i from 0 <= i < PS.degree:
     345        if PS.levels[i] <= PS.depth:
     346            ncells += 1
     347    return ncells
     348
     349cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end):
     350    """
     351    Makes sure that the first element of the segment of entries i with
     352    start <= i <= end is minimal.
     353    """
     354    cdef int i, min_loc = start, minimum = PS.entries[start]
     355    for i from start < i <= end:
     356        if PS.entries[i] < minimum:
     357            min_loc = i
     358            minimum = PS.entries[i]
     359    if min_loc != start:
     360        PS.entries[min_loc] = PS.entries[start]
     361        PS.entries[start] = minimum
     362
     363cdef inline bint PS_is_mcr(PartitionStack *PS, int m):
     364    """
     365    Returns whether PS.elements[m] (not m!) is the smallest element of its cell.
     366    """
     367    return m == 0 or PS.levels[m-1] <= PS.depth
     368
     369cdef inline bint PS_is_fixed(PartitionStack *PS, int m):
     370    """
     371    Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming
     372    PS_is_mcr(PS, m) is already true.
     373    """
     374    return PS.levels[m] <= PS.depth
     375
     376cdef inline int PS_clear(PartitionStack *PS):
     377    """
     378    Sets the current partition to the first shallower one, i.e. forgets about
     379    boundaries between cells that are new to the current level.
     380    """
     381    cdef int i, cur_start = 0
     382    for i from 0 <= i < PS.degree:
     383        if PS.levels[i] >= PS.depth:
     384            PS.levels[i] += 1
     385        if PS.levels[i] < PS.depth:
     386            PS_move_min_to_front(PS, cur_start, i)
     387            cur_start = i+1
     388
     389cdef inline int PS_move_all_mins_to_front(PartitionStack *PS):
     390    """
     391    Move minimal cell elements to the front of each cell.
     392    """
     393    cdef int i, cur_start = 0
     394    for i from 0 <= i < PS.degree:
     395        if PS.levels[i] <= PS.depth:
     396            PS_move_min_to_front(PS, cur_start, i)
     397            cur_start = i+1
     398
     399cdef inline int PS_split_point(PartitionStack *PS, int v):
     400    """
     401    Detaches the point v from the cell it is in, putting the singleton cell
     402    of just v in front. Returns the position where v is now located.
     403    """
     404    cdef int i = 0, index_of_v
     405    while PS.entries[i] != v:
     406        i += 1
     407    index_of_v = i
     408    while PS.levels[i] > PS.depth:
     409        i += 1
     410    if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \
     411       and PS.levels[index_of_v] > PS.depth:
     412        # if v is first (make sure v is not already alone)
     413        PS_move_min_to_front(PS, index_of_v+1, i)
     414        PS.levels[index_of_v] = PS.depth
     415        return index_of_v
     416    else:
     417        # If v is not at front, i.e. v is not minimal in its cell,
     418        # then move_min_to_front is not necessary since v will swap
     419        # with the first before making its own cell, leaving it at
     420        # the front of the other.
     421        i = index_of_v
     422        while i != 0 and PS.levels[i-1] > PS.depth:
     423            i -= 1
     424        PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v
     425        PS.entries[i+1] = PS.entries[i] # move the first (min) to second
     426        PS.entries[i] = v # place v first
     427        PS.levels[i] = PS.depth
     428        return i
     429
     430cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b):
     431    """
     432    Find the first occurrence of the smallest cell of size greater than one,
     433    store its entries to b, and return its minimum element.
     434    """
     435    cdef int i = 0, j = 0, location = 0, n = PS.degree
     436    bitset_zero(b)
     437    while 1:
     438        if PS.levels[i] <= PS.depth:
     439            if i != j and n > i - j + 1:
     440                n = i - j + 1
     441                location = j
     442            j = i + 1
     443        if PS.levels[i] == -1: break
     444        i += 1
     445    # location now points to the beginning of the first, smallest,
     446    # nontrivial cell
     447    i = location
     448    while 1:
     449        bitset_flip(b, PS.entries[i])
     450        if PS.levels[i] <= PS.depth: break
     451        i += 1
     452    return PS.entries[location]
     453
     454cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma):
     455    """
     456    Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i]
     457    denotes the entry of the ith cell of the discrete partition PS.
     458    """
     459    cdef int i
     460    for i from 0 <= i < PS1.degree:
     461        gamma[PS2.entries[i]] = PS1.entries[i]
     462
     463def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]):
     464    """
     465    Demonstration and testing.
     466   
     467    DOCTEST:
     468        sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \
     469        ... import PS_represent
     470        sage: PS_represent()
     471        Allocating PartitionStack...
     472        Allocation passed:
     473        (0,1,2,3,4,5,6,7,8,9)
     474        Checking that entries are in order and correct level.
     475        Everything seems in order, deallocating.
     476        Deallocated.
     477        Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]].
     478        PartitionStack's data:
     479        entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2]
     480        levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1]
     481        depth = 0, degree = 10
     482        (6|3,4,8,7|1,9,5|0,2)
     483        Checking PS_is_discrete:
     484        False
     485        Checking PS_num_cells:
     486        4
     487        Checking PS_is_mcr, min cell reps are:
     488        [6, 3, 1, 0]
     489        Checking PS_is_fixed, fixed elements are:
     490        [6]
     491        Copying PartitionStack:
     492        (6|3,4,8,7|1,9,5|0,2)
     493        Checking for consistency.
     494        Everything is consistent.
     495        Clearing copy:
     496        (0,3,4,8,7,1,9,5,6,2)
     497        Splitting point 6 from original:
     498        0
     499        (6|3,4,8,7|1,9,5|0,2)
     500        Splitting point 1 from original:
     501        5
     502        (6|3,4,8,7|1|5,9|0,2)
     503        Splitting point 8 from original:
     504        1
     505        (6|8|3,4,7|1|5,9|0,2)
     506        Splitting point 2 from original:
     507        8
     508        (6|8|3,4,7|1|5,9|2|0)
     509        Getting permutation from PS2->PS:
     510        [6, 1, 0, 8, 3, 9, 2, 7, 4, 5]
     511        Finding first smallest:
     512        Minimal element is 5, bitset is:
     513        0000010001
     514        Deallocating PartitionStacks.
     515        Done.
     516
     517    """
     518    cdef int i, n = sum([len(cell) for cell in partition]), *gamma
     519    cdef bitset_t b
     520    print "Allocating PartitionStack..."
     521    cdef PartitionStack *PS = PS_new(n, 1), *PS2
     522    if PS is NULL:
     523        print "Allocation failed!"
     524        return
     525    print "Allocation passed:"
     526    PS_print(PS)
     527    print "Checking that entries are in order and correct level."
     528    good = True
     529    for i from 0 <= i < n-1:
     530        if not (PS.entries[i] == i and PS.levels[i] == n):
     531            print "Failed at i = %d!"%i
     532            print PS.entries[i], PS.levels[i], i, n
     533            good = False
     534    if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1):
     535        print "Failed at i = %d!"%(n-1)
     536        good = False
     537    if not PS.degree == n or not PS.depth == 0:
     538        print "Incorrect degree or depth!"
     539        good = False
     540    if good: print "Everything seems in order, deallocating."
     541    PS_dealloc(PS)
     542    print "Deallocated."
     543    print "Creating PartitionStack from partition %s."%partition
     544    PS = PS_from_list(partition)
     545    print "PartitionStack's data:"
     546    print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n]
     547    print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n]
     548    print "depth = %d, degree = %d"%(PS.depth,PS.degree)
     549    PS_print(PS)
     550    print "Checking PS_is_discrete:"
     551    print "True" if PS_is_discrete(PS) else "False"
     552    print "Checking PS_num_cells:"
     553    print PS_num_cells(PS)
     554    print "Checking PS_is_mcr, min cell reps are:"
     555    L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)]
     556    print L
     557    print "Checking PS_is_fixed, fixed elements are:"
     558    print [PS.entries[l] for l in L if PS_is_fixed(PS, l)]
     559    print "Copying PartitionStack:"
     560    PS2 = PS_copy(PS)
     561    PS_print(PS2)
     562    print "Checking for consistency."
     563    good = True
     564    for i from 0 <= i < n:
     565        if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]:
     566            print "Failed at i = %d!"%i
     567            good = False
     568    if PS.degree != PS2.degree or PS.depth != PS2.depth:
     569        print "Failure with degree or depth!"
     570        good = False
     571    if good:
     572        print "Everything is consistent."
     573    print "Clearing copy:"
     574    PS_clear(PS2)
     575    PS_print(PS2)
     576    for s in splits:
     577        print "Splitting point %d from original:"%s
     578        print PS_split_point(PS, s)
     579        PS_print(PS)
     580    print "Getting permutation from PS2->PS:"
     581    gamma = <int *> sage_malloc(n * sizeof(int))
     582    PS_get_perm_from(PS, PS2, gamma)
     583    print [gamma[i] for i from 0 <= i < n]
     584    sage_free(gamma)
     585    print "Finding first smallest:"
     586    bitset_init(b, n)
     587    i = PS_first_smallest(PS, b)
     588    print "Minimal element is %d, bitset is:"%i
     589    print bitset_string(b)
     590    bitset_clear(b)
     591    print "Deallocating PartitionStacks."
     592    PS_dealloc(PS)
     593    PS_dealloc(PS2)
     594    print "Done."
     595
     596# Functions
     597
     598cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S,
     599    int (*refine_and_return_invariant)\
     600         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
     601    int *cells_to_refine_by):
     602    """
     603    Make the partition stack one longer by copying the last partition in the
     604    stack, split off a given point, and refine. Return the invariant given by
     605    the refinement function.
     606   
     607    INPUT:
     608    PS -- the partition stack to refine
     609    v -- the point to split
     610    S -- the structure
     611    refine_and_return_invariant -- the refinement function provided
     612    cells_to_refine_by -- an array, contents ignored
     613
     614    """
     615    PS.depth += 1
     616    PS_clear(PS)
     617    cells_to_refine_by[0] = PS_split_point(PS, v)
     618    return refine_and_return_invariant(PS, S, cells_to_refine_by, 1)
     619
  • new file sage/groups/perm_gps/partn_ref/double_coset.pxd

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/double_coset.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 'data_structures_pxd.pxi' # includes bitsets
     12
     13from sage.rings.integer cimport Integer
     14
     15cdef int *double_coset( object, object, int **, int *, int,
     16    bint (*)(PartitionStack *, object),
     17    int (*)(PartitionStack *, object, int *, int),
     18    int (*)(int *, int *, object, object))
     19
     20
     21
  • new file sage/groups/perm_gps/partn_ref/double_coset.pyx

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/double_coset.pyx
    - +  
     1r"""
     2Double cosets
     3
     4This module implements a general algorithm for computing double coset problems
     5for pairs of objects. The class of objects in question must be some kind
     6of structure for which an isomorphism is a permutation in $S_n$ for some $n$,
     7which we call here the order of the object. Given objects $X$ and $Y$,
     8the program returns an isomorphism in list permutation form if $X \cong Y$, and
     9a NULL pointer otherwise.
     10
     11In order to take advantage of the algorithms in this module for a specific kind
     12of object, one must implement (in Cython) three functions which will be specific
     13to the kind of objects in question. Pointers to these functions are passed to
     14the main function of the module, which is \code{double_coset}. For specific
     15examples of implementations of these functions, see any of the files in
     16\code{sage.groups.perm_gps.partn_ref} beginning with "refinement." They are:
     17
     18A. \code{refine_and_return_invariant}:
     19   
     20    Signature:
     21   
     22    \code{int refine_and_return_invariant(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len)}
     23   
     24    This function should split up cells in the partition at the top of the
     25    partition stack in such a way that any automorphism that respects the
     26    partition also respects the resulting partition. The array
     27    cells_to_refine_by is a list of the beginning positions of some cells which
     28    have been changed since the last refinement. It is not necessary to use
     29    this in an implementation of this function, but it will affect performance.
     30    One should consult \code{refinement_graphs} for more details and ideas for
     31    particular implementations.
     32   
     33    Output:
     34   
     35    An integer $I$ invariant under the orbits of $S_n$.  That is, if
     36    $\gamma \in S_n$, then
     37    $$ I(G, PS, cells_to_refine_by) = I( \gamma(G), \gamma(PS), \gamma(cells_to_refine_by) ) .$$
     38
     39
     40B. \code{compare_structures}:
     41   
     42    Signature:
     43   
     44    \code{int compare_structures(int *gamma_1, int *gamma_2, object S)}
     45   
     46    This function must implement a total ordering on the set of objects of fixed
     47    order. Return -1 if \code{gamma_1(S) < gamma_2(S)}, 0 if
     48    \code{gamma_1(S) == gamma_2(S)}, 1 if \code{gamma_1(S) > gamma_2(S)}.
     49
     50C. \code{all_children_are_equivalent}:
     51   
     52    Signature:
     53   
     54    \code{bint all_children_are_equivalent(PartitionStack *PS, object S)}
     55   
     56    This function must return False unless it is the case that each discrete
     57    partition finer than the top of the partition stack is equivalent to the
     58    others under some automorphism of S. The converse need not hold: if this is
     59    indeed the case, it still may return False. This function is originally used
     60    as a consequence of Lemma 2.25 in [1].
     61
     62DOCTEST:
     63    sage: import sage.groups.perm_gps.partn_ref.double_coset
     64
     65REFERENCE:
     66
     67    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
     68        Vol. 30 (1981), pp. 45-87.
     69
     70    [2] Leon, Jeffrey. Permutation Group Algorithms Based on Partitions, I:
     71        Theory and Algorithms. J. Symbolic Computation, Vol. 12 (1991), pp.
     72        533-583.
     73
     74"""
     75
     76#*****************************************************************************
     77#      Copyright (C) 2006 - 2008 Robert L. Miller <rlmillster@gmail.com>
     78#
     79# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
     80#                         http://www.gnu.org/licenses/
     81#*****************************************************************************
     82
     83include 'data_structures_pyx.pxi' # includes bitsets
     84
     85cdef inline int cmp(int a, int b):
     86    if a < b: return -1
     87    elif a == b: return 0
     88    else: return 1
     89
     90cdef inline bint stacks_are_equivalent(PartitionStack *PS1, PartitionStack *PS2):
     91    cdef int i, j, depth = min(PS1.depth, PS2.depth)
     92    for i from 0 <= i < PS1.degree:
     93        j = cmp(PS1.levels[i], PS2.levels[i])
     94        if j == 0: continue
     95        if ( (j < 0 and PS1.levels[i] <= depth and PS2.levels[i] > depth)
     96            or (j > 0 and PS2.levels[i] <= depth and PS1.levels[i] > depth) ):
     97            return 0
     98    return 1
     99
     100# Functions
     101
     102cdef int *double_coset(object S1, object S2, int **partition1, int *ordering2,
     103    int n, bint (*all_children_are_equivalent)(PartitionStack *PS, object S),
     104    int (*refine_and_return_invariant)\
     105         (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len),
     106    int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2) ):
     107    """
     108    Traverse the search space for double coset calculation.
     109   
     110    INPUT:
     111    S1, S2 -- pointers to the structures
     112    partition1 -- array representing a partition of the points of S1
     113    ordering2 -- an ordering of the points of S2 representing a second partition
     114    n -- the number of points (points are assumed to be 0,1,...,n-1)
     115    all_children_are_equivalent -- pointer to a function
     116        INPUT:
     117        PS -- pointer to a partition stack
     118        S -- pointer to the structure
     119        OUTPUT:
     120        bint -- returns True if it can be determined that all refinements below
     121            the current one will result in an equivalent discrete partition
     122    refine_and_return_invariant -- pointer to a function
     123        INPUT:
     124        PS -- pointer to a partition stack
     125        S -- pointer to the structure
     126        alpha -- an array consisting of numbers, which indicate the starting
     127            positions of the cells to refine against (will likely be modified)
     128        OUTPUT:
     129        int -- returns an invariant under application of arbitrary permutations
     130    compare_structures -- pointer to a function
     131        INPUT:
     132        gamma_1, gamma_2 -- (list) permutations of the points of S1 and S2
     133        S1, S2 -- pointers to the structures
     134        OUTPUT:
     135        int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp),
     136            such that the set of all structures is well-ordered
     137    OUTPUT:
     138    If S1 and S2 are isomorphic, a pointer to an integer array representing an
     139    isomorphism. Otherwise, a NULL pointer.
     140
     141    """
     142    cdef PartitionStack *current_ps, *first_ps, *left_ps
     143    cdef int first_meets_current = -1
     144    cdef int current_kids_are_same = 1
     145    cdef int first_kids_are_same
     146   
     147    cdef int *indicators
     148   
     149    cdef OrbitPartition *orbits_of_subgroup
     150    cdef int subgroup_primary_orbit_size = 0
     151    cdef int minimal_in_primary_orbit
     152   
     153    cdef bitset_t *fixed_points_of_generators # i.e. fp
     154    cdef bitset_t *minimal_cell_reps_of_generators # i.e. mcr
     155    cdef int len_of_fp_and_mcr = 100
     156    cdef int index_in_fp_and_mcr = -1
     157   
     158    cdef bitset_t *vertices_to_split
     159    cdef bitset_t vertices_have_been_reduced
     160    cdef int *permutation, *id_perm, *cells_to_refine_by
     161    cdef int *vertices_determining_current_stack
     162   
     163    cdef int *output
     164   
     165    cdef int i, j, k
     166    cdef bint discrete, automorphism, update_label
     167    cdef bint backtrack, new_vertex, narrow, mem_err = 0
     168   
     169    if n == 0:
     170        return NULL
     171   
     172    indicators = <int *> sage_malloc(n * sizeof(int))
     173   
     174    fixed_points_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
     175    minimal_cell_reps_of_generators = <bitset_t *> sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) )
     176   
     177    vertices_to_split = <bitset_t *> sage_malloc( n * sizeof(bitset_t) )
     178    permutation = <int *> sage_malloc( n * sizeof(int) )
     179    id_perm = <int *> sage_malloc( n * sizeof(int) )
     180    cells_to_refine_by = <int *> sage_malloc( n * sizeof(int) )
     181    vertices_determining_current_stack = <int *> sage_malloc( n * sizeof(int) )
     182
     183    output = <int *> sage_malloc( n * sizeof(int) )
     184   
     185    current_ps = PS_new(n, 0)
     186    left_ps = PS_new(n, 0)
     187    orbits_of_subgroup = OP_new(n)
     188   
     189    # Check for allocation failures:
     190    if indicators                         is NULL or \
     191       fixed_points_of_generators         is NULL or \
     192       minimal_cell_reps_of_generators    is NULL or \
     193       vertices_to_split                  is NULL or \
     194       permutation                        is NULL or \
     195       id_perm                            is NULL or \
     196       cells_to_refine_by                 is NULL or \
     197       vertices_determining_current_stack is NULL or \
     198       current_ps                         is NULL or \
     199       orbits_of_subgroup                 is NULL or \
     200       output                             is NULL:
     201        if indicators                         is not NULL:
     202            sage_free(indicators)
     203        if fixed_points_of_generators         is not NULL:
     204            sage_free(fixed_points_of_generators)
     205        if minimal_cell_reps_of_generators    is not NULL:
     206            sage_free(minimal_cell_reps_of_generators)
     207        if vertices_to_split                  is not NULL:
     208            sage_free(vertices_to_split)
     209        if permutation                        is not NULL:
     210            sage_free(permutation)
     211        if id_perm                            is not NULL:
     212            sage_free(id_perm)
     213        if cells_to_refine_by                 is not NULL:
     214            sage_free(cells_to_refine_by)
     215        if vertices_determining_current_stack is not NULL:
     216            sage_free(vertices_determining_current_stack)
     217        if current_ps is not NULL:
     218            PS_dealloc(current_ps)
     219        if orbits_of_subgroup is not NULL:
     220            OP_dealloc(orbits_of_subgroup)
     221        if output is not NULL:
     222            sage_free(output)
     223        raise MemoryError
     224   
     225    # Initialize bitsets, checking for allocation failures:
     226    cdef bint succeeded = 1
     227    for i from 0 <= i < len_of_fp_and_mcr:
     228        try:
     229            bitset_init(fixed_points_of_generators[i], n)
     230        except MemoryError:
     231            succeeded = 0
     232            for j from 0 <= j < i:
     233                bitset_clear(fixed_points_of_generators[j])
     234                bitset_clear(minimal_cell_reps_of_generators[j])
     235            break
     236        try:
     237            bitset_init(minimal_cell_reps_of_generators[i], n)
     238        except MemoryError:
     239            succeeded = 0
     240            for j from 0 <= j < i:
     241                bitset_clear(fixed_points_of_generators[j])
     242                bitset_clear(minimal_cell_reps_of_generators[j])
     243            bitset_clear(fixed_points_of_generators[i])
     244            break
     245    if succeeded:
     246        for i from 0 <= i < n:
     247            try:
     248                bitset_init(vertices_to_split[i], n)
     249            except MemoryError:
     250                succeeded = 0
     251                for j from 0 <= j < i:
     252                    bitset_clear(vertices_to_split[j])
     253                for j from 0 <= j < len_of_fp_and_mcr:
     254                    bitset_clear(fixed_points_of_generators[j])
     255                    bitset_clear(minimal_cell_reps_of_generators[j])
     256                break
     257    if succeeded:
     258        try:
     259            bitset_init(vertices_have_been_reduced, n)
     260        except MemoryError:
     261            succeeded = 0
     262            for j from 0 <= j < n:
     263                bitset_clear(vertices_to_split[j])
     264            for j from 0 <= j < len_of_fp_and_mcr:
     265                bitset_clear(fixed_points_of_generators[j])
     266                bitset_clear(minimal_cell_reps_of_generators[j])
     267    if not succeeded:
     268        sage_free(indicators)
     269        sage_free(fixed_points_of_generators)
     270        sage_free(minimal_cell_reps_of_generators)
     271        sage_free(vertices_to_split)
     272        sage_free(permutation)
     273        sage_free(id_perm)
     274        sage_free(cells_to_refine_by)
     275        sage_free(vertices_determining_current_stack)
     276        PS_dealloc(current_ps)
     277        PS_dealloc(left_ps)
     278        OP_dealloc(orbits_of_subgroup)
     279        raise MemoryError
     280   
     281    bitset_zero(vertices_have_been_reduced)
     282   
     283    # set up the identity permutation
     284    for i from 0 <= i < n:
     285        id_perm[i] = i
     286    if ordering2 is NULL:
     287        ordering2 = id_perm
     288   
     289    # Copy data of partition to left_ps, and
     290    # ordering of that partition to current_ps.
     291    i = 0
     292    j = 0
     293    while partition1[i] is not NULL:
     294        k = 0
     295        while partition1[i][k] != -1:
     296            left_ps.entries[j+k] = partition1[i][k]
     297            left_ps.levels[j+k] = n
     298            current_ps.entries[j+k] = ordering2[j+k]
     299            current_ps.levels[j+k] = n
     300            k += 1
     301        left_ps.levels[j+k-1] = 0
     302        current_ps.levels[j+k-1] = 0
     303        PS_move_min_to_front(current_ps, j, j+k-1)
     304        j += k
     305        i += 1
     306    left_ps.levels[j-1] = -1
     307    current_ps.levels[j-1] = -1
     308    left_ps.depth = 0
     309    current_ps.depth = 0
     310    left_ps.degree = n
     311    current_ps.degree = n
     312   
     313    # default values of "infinity"
     314    for i from 0 <= i < n:
     315        indicators[i] = -1
     316   
     317    cdef bint possible = 1
     318    cdef bint unknown = 1
     319   
     320    # Our first refinement needs to check every cell of the partition,
     321    # so cells_to_refine_by needs to be a list of pointers to each cell.
     322    j = 1
     323    cells_to_refine_by[0] = 0
     324    for i from 0 < i < n:
     325        if left_ps.levels[i-1] == 0:
     326            cells_to_refine_by[j] = i
     327            j += 1
     328   
     329    k = refine_and_return_invariant(left_ps, S1, cells_to_refine_by, j)
     330   
     331    j = 1
     332    cells_to_refine_by[0] = 0
     333    for i from 0 < i < n:
     334        if current_ps.levels[i-1] == 0:
     335            cells_to_refine_by[j] = i
     336            j += 1
     337    j = refine_and_return_invariant(current_ps, S2, cells_to_refine_by, j)
     338    if k != j:
     339        possible = 0; unknown = 0
     340    elif not stacks_are_equivalent(left_ps, current_ps):
     341        possible = 0; unknown = 0
     342    else:
     343        PS_move_all_mins_to_front(current_ps)
     344
     345    first_ps = NULL
     346    # Refine down to a discrete partition
     347    while not PS_is_discrete(left_ps) and possible:
     348        k = PS_first_smallest(left_ps, vertices_to_split[left_ps.depth])
     349        i = left_ps.depth
     350        indicators[i] = split_point_and_refine(left_ps, k, S1, refine_and_return_invariant, cells_to_refine_by)
     351        vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
     352        bitset_unset(vertices_have_been_reduced, current_ps.depth)
     353        while 1:
     354            j =  split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by)
     355            if indicators[i] != j:
     356                possible = 0
     357            elif not stacks_are_equivalent(left_ps, current_ps):
     358                possible = 0
     359            else:
     360                PS_move_all_mins_to_front(current_ps)
     361            if not possible:
     362                j = vertices_determining_current_stack[i] + 1
     363                j = bitset_next(vertices_to_split[i], j)
     364                if j == -1:
     365                    break
     366                else:
     367                    possible = 1
     368                    vertices_determining_current_stack[i] = j
     369                    current_ps.depth -= 1 # reset for next refinement
     370            else: break
     371        if not all_children_are_equivalent(current_ps, S2):
     372            current_kids_are_same = current_ps.depth + 1
     373   
     374    if possible:
     375        if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
     376            unknown = 0
     377        else:
     378            first_meets_current = current_ps.depth
     379            first_kids_are_same = current_ps.depth
     380            first_ps = PS_copy(current_ps)
     381            current_ps.depth -= 1
     382   
     383    # Main loop:
     384    while possible and unknown and current_ps.depth != -1:
     385       
     386        # I. Search for a new vertex to split, and update subgroup information
     387        new_vertex = 0
     388        if current_ps.depth > first_meets_current:
     389            # If we are not at a node of the first stack, reduce size of
     390            # vertices_to_split by using the symmetries we already know.
     391            if not bitset_check(vertices_have_been_reduced, current_ps.depth):
     392                for i from 0 <= i <= index_in_fp_and_mcr:
     393                    j = 0
     394                    while j < current_ps.depth and bitset_check(fixed_points_of_generators[i], vertices_determining_current_stack[j]):
     395                        j += 1
     396                    # If each vertex split so far is fixed by generator i,
     397                    # then remove elements of vertices_to_split which are
     398                    # not minimal in their orbits under generator i.
     399                    if j == current_ps.depth:
     400                        for k from 0 <= k < n:
     401                            if bitset_check(vertices_to_split[current_ps.depth], k) and not bitset_check(minimal_cell_reps_of_generators[i], k):
     402                                bitset_flip(vertices_to_split[current_ps.depth], k)
     403                bitset_flip(vertices_have_been_reduced, current_ps.depth)
     404            # Look for a new point to split.
     405            i = vertices_determining_current_stack[current_ps.depth] + 1
     406            i = bitset_next(vertices_to_split[current_ps.depth], i)
     407            if i != -1:
     408                # There is a new point.
     409                vertices_determining_current_stack[current_ps.depth] = i
     410                new_vertex = 1
     411            else:
     412                # No new point: backtrack.
     413                current_ps.depth -= 1
     414        else:
     415            # If we are at a node of the first stack, the above reduction
     416            # will not help. Also, we must update information about
     417            # primary orbits here.
     418            if current_ps.depth < first_meets_current:
     419                # If we are done searching under this part of the first
     420                # stack, then first_meets_current is one higher, and we
     421                # are looking at a new primary orbit (corresponding to a
     422                # larger subgroup in the stabilizer chain).
     423                first_meets_current = current_ps.depth
     424                for i from 0 <= i < n:
     425                    if bitset_check(vertices_to_split[current_ps.depth], i):
     426                        minimal_in_primary_orbit = i
     427                        break
     428            while 1:
     429                i = vertices_determining_current_stack[current_ps.depth]
     430                # This was the last point to be split here.
     431                # If it is in the same orbit as minimal_in_primary_orbit,
     432                # then count it as an element of the primary orbit.
     433                if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit):
     434                    subgroup_primary_orbit_size += 1
     435                # Look for a new point to split.
     436                i += 1
     437                i = bitset_next(vertices_to_split[current_ps.depth], i)
     438                if i != -1:
     439                    # There is a new point.
     440                    vertices_determining_current_stack[current_ps.depth] = i
     441                    if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i:
     442                        new_vertex = 1
     443                        break
     444                else:
     445                    # No new point: backtrack.
     446                    # Note that now, we are backtracking up the first stack.
     447                    vertices_determining_current_stack[current_ps.depth] = -1
     448                    # If every choice of point to split gave something in the
     449                    # (same!) primary orbit, then all children of the first
     450                    # stack at this point are equivalent.
     451                    j = 0
     452                    for i from 0 <= i < n:
     453                        if bitset_check(vertices_to_split[current_ps.depth], i):
     454                            j += 1
     455                    if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1:
     456                        first_kids_are_same = current_ps.depth
     457                    # Backtrack.
     458                    subgroup_primary_orbit_size = 0
     459                    current_ps.depth -= 1
     460                    break
     461        if not new_vertex:
     462            continue
     463   
     464        if current_kids_are_same > current_ps.depth + 1:
     465            current_kids_are_same = current_ps.depth + 1
     466       
     467        # II. Refine down to a discrete partition, or until
     468        # we leave the part of the tree we are interested in
     469        discrete = 0
     470        while 1:
     471            i = current_ps.depth
     472            while 1:
     473                k = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by)
     474                PS_move_all_mins_to_front(current_ps)
     475                if indicators[i] != k:
     476                    possible = 0
     477                elif not stacks_are_equivalent(left_ps, current_ps):
     478                    possible = 0
     479                if not possible:
     480                    j = vertices_determining_current_stack[i] + 1
     481                    j = bitset_next(vertices_to_split[i], j)
     482                    if j == -1:
     483                        break
     484                    else:
     485                        possible = 1
     486                        vertices_determining_current_stack[i] = j
     487                        current_ps.depth -= 1 # reset for next refinement
     488                else: break
     489            if not possible:
     490                break
     491            if PS_is_discrete(current_ps):
     492                break
     493            vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth])
     494            bitset_unset(vertices_have_been_reduced, current_ps.depth)
     495            if not all_children_are_equivalent(current_ps, S2):
     496                current_kids_are_same = current_ps.depth + 1
     497       
     498        # III. Check for automorphisms and isomorphisms
     499        automorphism = 0
     500        if possible:
     501            PS_get_perm_from(current_ps, first_ps, permutation)
     502            if compare_structures(permutation, id_perm, S2, S2) == 0:
     503                automorphism = 1
     504        if not automorphism and possible:
     505            # if we get here, discrete must be true
     506            if current_ps.depth != left_ps.depth:
     507                possible = 0
     508            elif compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0:
     509                unknown = 0
     510                break # main loop
     511            else:
     512                possible = 0
     513        if automorphism:
     514            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
     515                index_in_fp_and_mcr += 1
     516            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
     517            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     518            for i from 0 <= i < n:
     519                if permutation[i] == i:
     520                    bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
     521                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     522                else:
     523                    bitset_unset(fixed_points_of_generators[index_in_fp_and_mcr], i)
     524                    k = i
     525                    j = permutation[i]
     526                    while j != i:
     527                        if j < k: k = j
     528                        j = permutation[j]
     529                    if k == i:
     530                        bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     531                    else:
     532                        bitset_unset(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     533            current_ps.depth = first_meets_current
     534            if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser
     535                if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit:
     536                    continue # main loop
     537            if bitset_check(vertices_have_been_reduced, current_ps.depth):
     538                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     539            continue # main loop
     540        if not possible:
     541            possible = 1
     542            i = current_ps.depth
     543            current_ps.depth = min(first_kids_are_same-1, current_kids_are_same-1)
     544            if i == current_kids_are_same:
     545                continue # main loop
     546            if index_in_fp_and_mcr < len_of_fp_and_mcr - 1:
     547                index_in_fp_and_mcr += 1
     548            bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr])
     549            bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     550            j = current_ps.depth
     551            current_ps.depth = i # just for mcr and fixed functions...
     552            for i from 0 <= i < n:
     553                if PS_is_mcr(current_ps, i):
     554                    bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i)
     555                    if PS_is_fixed(current_ps, i):
     556                        bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i)
     557            current_ps.depth = j
     558            if bitset_check(vertices_have_been_reduced, current_ps.depth):
     559                bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr])
     560   
     561    # End of main loop.
     562   
     563    if possible and not unknown:
     564        for i from 0 <= i < n:
     565            output[left_ps.entries[i]] = current_ps.entries[i]
     566    else:
     567        sage_free(output)
     568        output = NULL
     569   
     570    # Deallocate:
     571    for i from 0 <= i < len_of_fp_and_mcr:
     572        bitset_clear(fixed_points_of_generators[i])
     573        bitset_clear(minimal_cell_reps_of_generators[i])
     574    for i from 0 <= i < n:
     575        bitset_clear(vertices_to_split[i])
     576    bitset_clear(vertices_have_been_reduced)
     577    sage_free(indicators)
     578    sage_free(fixed_points_of_generators)
     579    sage_free(minimal_cell_reps_of_generators)
     580    sage_free(vertices_to_split)
     581    sage_free(permutation)
     582    sage_free(id_perm)
     583    sage_free(cells_to_refine_by)
     584    sage_free(vertices_determining_current_stack)
     585    PS_dealloc(current_ps)
     586    PS_dealloc(left_ps)
     587    OP_dealloc(orbits_of_subgroup)
     588    if first_ps is not NULL: PS_dealloc(first_ps)
     589   
     590    return output
     591
     592
     593
     594
     595
  • sage/groups/perm_gps/partn_ref/refinement_binary.pxd

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_binary.pxd
    a b  
    88
    99include '../../../ext/cdefs.pxi'
    1010include '../../../ext/stdsage.pxi'
    11 include '../../../misc/bitset_pxd.pxi'
     11include 'data_structures_pxd.pxi' # includes bitsets
    1212
    1313from sage.rings.integer cimport Integer
    14 from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport PartitionStack, PS_print, PS_new, PS_dealloc, PS_clear, PS_is_discrete, PS_move_min_to_front, PS_num_cells, get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     14from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     15from double_coset cimport double_coset
    1516
    1617cdef class BinaryCodeStruct:
    1718    cdef bitset_s *alpha_is_wd # single bitset of length nwords + degree
     
    3637cdef int ith_word_nonlinear(BinaryCodeStruct, int, bitset_s *)
    3738
    3839cdef int refine_by_bip_degree(PartitionStack *, object, int *, int)
    39 cdef int compare_linear_codes(int *, int *, object)
    40 cdef int compare_nonlinear_codes(int *, int *, object)
     40cdef int compare_linear_codes(int *, int *, object, object)
     41cdef int compare_nonlinear_codes(int *, int *, object, object)
    4142cdef bint all_children_are_equivalent(PartitionStack *, object)
    4243cdef inline int word_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *)
    4344cdef inline int col_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *)
  • sage/groups/perm_gps/partn_ref/refinement_binary.pyx

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_binary.pyx
    a b  
    2222#                         http://www.gnu.org/licenses/
    2323#*****************************************************************************
    2424
    25 include '../../../misc/bitset.pxi'
     25include 'data_structures_pyx.pxi' # includes bitsets
     26
    2627from sage.matrix.matrix import is_Matrix
    2728
    2829cdef class LinearBinaryCodeStruct(BinaryCodeStruct):
     
    9596       
    9697        self.output = NULL
    9798        self.ith_word = &ith_word_linear
    98         self.first_time = 1
    9999   
    100100    def run(self, partition=None):
    101101        """
     
    238238                part[i][j] = partition[i][j]
    239239            part[i][len(partition[i])] = -1
    240240        part[len(partition)] = NULL
     241        self.first_time = 1
    241242       
    242243        self.output = get_aut_gp_and_can_lab(self, part, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, 1, 1, 1)
    243244       
     
    295296        if self.output is NULL:
    296297            self.run()
    297298        return [self.output.relabeling[i] for i from 0 <= i < self.degree]
     299   
     300    def is_isomorphic(self, LinearBinaryCodeStruct other):
     301        """
     302        Calculate whether self is isomorphic to other.
     303       
     304        EXAMPLES:
     305            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct
     306
     307            sage: B = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
     308            sage: C = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,1],[1,1,0,1,1,0]]))
     309            sage: B.is_isomorphic(C)
     310            [0, 1, 2, 5, 3, 4]
     311           
     312        """
     313        cdef int **part, i, j
     314        cdef int *output, *ordering
     315        partition = [range(self.degree)]
     316        part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     317        ordering = <int *> sage_malloc(self.degree * sizeof(int))
     318        if part is NULL or ordering is NULL:
     319            if part is not NULL: sage_free(part)
     320            if ordering is not NULL: sage_free(ordering)
     321            raise MemoryError
     322        for i from 0 <= i < len(partition):
     323            part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
     324            if part[i] is NULL:
     325                for j from 0 <= j < i:
     326                    sage_free(part[j])
     327                sage_free(part)
     328                raise MemoryError
     329            for j from 0 <= j < len(partition[i]):
     330                part[i][j] = partition[i][j]
     331            part[i][len(partition[i])] = -1
     332        part[len(partition)] = NULL
     333        for i from 0 <= i < self.degree:
     334            ordering[i] = i
     335        self.first_time = 1
     336        other.first_time = 1
     337       
     338        output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes)
     339       
     340        for i from 0 <= i < len(partition):
     341            sage_free(part[i])
     342        sage_free(part)
     343        sage_free(ordering)
     344
     345        if output is NULL:
     346            return False
     347        else:
     348            output_py = [output[i] for i from 0 <= i < self.degree]
     349            sage_free(output)
     350            return output_py
    298351   
    299352    def __dealloc__(self):
    300353        cdef int j
     
    402455       
    403456        self.output = NULL
    404457        self.ith_word = &ith_word_nonlinear
    405         self.first_time = 1
    406458
    407459    def __dealloc__(self):
    408460        cdef int j
     
    477529                part[i][j] = partition[i][j]
    478530            part[i][len(partition[i])] = -1
    479531        part[len(partition)] = NULL
     532        self.first_time = 1
    480533       
    481534        self.output = get_aut_gp_and_can_lab(self, part, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, 1, 1, 1)
    482535       
     
    535588            self.run()
    536589        return [self.output.relabeling[i] for i from 0 <= i < self.degree]
    537590
     591    def is_isomorphic(self, NonlinearBinaryCodeStruct other):
     592        """
     593        Calculate whether self is isomorphic to other.
     594       
     595        EXAMPLES:
     596            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
     597
     598            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
     599            sage: C = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,0,0,1,1],[1,1,1,1,0,0]]))
     600            sage: B.is_isomorphic(C)
     601            [2, 3, 0, 1, 4, 5]
     602           
     603        """
     604        cdef int **part, i, j
     605        cdef int *output, *ordering
     606        partition = [range(self.degree)]
     607        part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     608        ordering = <int *> sage_malloc(self.degree * sizeof(int))
     609        if part is NULL or ordering is NULL:
     610            if part is not NULL: sage_free(part)
     611            if ordering is not NULL: sage_free(ordering)
     612            raise MemoryError
     613        for i from 0 <= i < len(partition):
     614            part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
     615            if part[i] is NULL:
     616                for j from 0 <= j < i:
     617                    sage_free(part[j])
     618                sage_free(part)
     619                raise MemoryError
     620            for j from 0 <= j < len(partition[i]):
     621                part[i][j] = partition[i][j]
     622            part[i][len(partition[i])] = -1
     623        part[len(partition)] = NULL
     624        for i from 0 <= i < self.degree:
     625            ordering[i] = i
     626        self.first_time = 1
     627        other.first_time = 1
     628       
     629        output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes)
     630       
     631        for i from 0 <= i < len(partition):
     632            sage_free(part[i])
     633        sage_free(part)
     634        sage_free(ordering)
     635
     636        if output is NULL:
     637            return False
     638        else:
     639            output_py = [output[i] for i from 0 <= i < self.degree]
     640            sage_free(output)
     641            return output_py
     642   
    538643cdef int ith_word_nonlinear(BinaryCodeStruct self, int i, bitset_s *word):
    539644    cdef NonlinearBinaryCodeStruct NBCS = <NonlinearBinaryCodeStruct> self
    540645    bitset_copy(word, &NBCS.words[i])
     
    606711                if necessary_to_split_cell:
    607712                    invariant += 8
    608713                    first_largest_subcell = sort_by_function(col_ps, current_cell, col_degrees, col_counts, col_output, BCS.nwords+1)
     714                    invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps)
    609715                    invariant += first_largest_subcell
    610716                    against_index = current_cell_against
    611717                    while against_index < ctrb_len:
     
    622728                        r += 1
    623729                        if r >= i:
    624730                            break
    625                     invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps)
    626731                    invariant += (i - current_cell)
    627732                current_cell = i
    628733        else:
     
    641746                if necessary_to_split_cell:
    642747                    invariant += 64
    643748                    first_largest_subcell = sort_by_function(word_ps, current_cell, word_degrees, word_counts, word_output, BCS.degree+1)
     749                    invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps)
    644750                    invariant += first_largest_subcell
    645751                    against_index = current_cell_against
    646752                    while against_index < ctrb_len:
     
    658764                        r += 1
    659765                        if r >= i:
    660766                            break
    661                     invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps)
    662767                    invariant += (i - current_cell)
    663768                current_cell = i
    664769        current_cell_against += 1
    665770    return invariant
    666771
    667 cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S):
     772cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S1, object S2):
    668773    """
    669     Compare gamma_1(B) and gamma_2(B).
     774    Compare gamma_1(S1) and gamma_2(S2).
    670775   
    671     Return return -1 if gamma_1(B) < gamma_2(B), 0 if gamma_1(B) == gamma_2(B),
    672     1 if gamma_1(B) > gamma_2(B).  (Just like the python \code{cmp}) function.
     776    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
     777    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
    673778   
    674779    Abstractly, what this function does is relabel the basis of B by gamma_1 and
    675780    gamma_2, run a row reduction on each, and verify that the matrices are the
     
    680785   
    681786    INPUT:
    682787    gamma_1, gamma_2 -- list permutations (inverse)
    683     S -- a binary code struct object
     788    S1, S2 -- binary code struct objects
    684789   
    685790    """
    686791    cdef int i, piv_loc_1, piv_loc_2, cur_col, cur_row=0
    687792    cdef bint is_pivot_1, is_pivot_2
    688     cdef LinearBinaryCodeStruct BCS = <LinearBinaryCodeStruct> S
    689     cdef bitset_s *basis_1 = BCS.scratch_bitsets                   # len = dim
    690     cdef bitset_s *basis_2 = &BCS.scratch_bitsets[BCS.dimension]   # len = dim
    691     cdef bitset_s *pivots  = &BCS.scratch_bitsets[2*BCS.dimension] # len 1
    692     cdef bitset_s *temp  = &BCS.scratch_bitsets[2*BCS.dimension+1] # len 1
    693     for i from 0 <= i < BCS.dimension:
    694         bitset_copy(&basis_1[i], &BCS.basis[i])
    695         bitset_copy(&basis_2[i], &BCS.basis[i])
     793    cdef LinearBinaryCodeStruct BCS1 = <LinearBinaryCodeStruct> S1
     794    cdef LinearBinaryCodeStruct BCS2 = <LinearBinaryCodeStruct> S2
     795    cdef bitset_s *basis_1 = BCS1.scratch_bitsets                   # len = dim
     796    cdef bitset_s *basis_2 = &BCS1.scratch_bitsets[BCS1.dimension]   # len = dim
     797    cdef bitset_s *pivots  = &BCS1.scratch_bitsets[2*BCS1.dimension] # len 1
     798    cdef bitset_s *temp  = &BCS1.scratch_bitsets[2*BCS1.dimension+1] # len 1
     799    for i from 0 <= i < BCS1.dimension:
     800        bitset_copy(&basis_1[i], &BCS1.basis[i])
     801        bitset_copy(&basis_2[i], &BCS2.basis[i])
    696802    bitset_zero(pivots)
    697     for cur_col from 0 <= cur_col < BCS.degree:
     803    for cur_col from 0 <= cur_col < BCS1.degree:
    698804        is_pivot_1 = 0
    699805        is_pivot_2 = 0
    700         for i from cur_row <= i < BCS.dimension:
     806        for i from cur_row <= i < BCS1.dimension:
    701807            if bitset_check(&basis_1[i], gamma_1[cur_col]):
    702808                is_pivot_1 = 1
    703809                piv_loc_1 = i
    704810                break
    705         for i from cur_row <= i < BCS.dimension:
     811        for i from cur_row <= i < BCS1.dimension:
    706812            if bitset_check(&basis_2[i], gamma_2[cur_col]):
    707813                is_pivot_2 = 1
    708814                piv_loc_2 = i
     
    724830                    bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
    725831                if bitset_check(&basis_2[i], gamma_2[cur_col]):
    726832                    bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row])
    727             for i from cur_row < i < BCS.dimension:
     833            for i from cur_row < i < BCS1.dimension:
    728834                if bitset_check(&basis_1[i], gamma_1[cur_col]):
    729835                    bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
    730836                if bitset_check(&basis_2[i], gamma_2[cur_col]):
     
    736842                    return <int>bitset_check(&basis_2[i], gamma_2[cur_col]) - <int>bitset_check(&basis_1[i], gamma_1[cur_col])
    737843    return 0
    738844
    739 cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S):
     845cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S1, object S2):
    740846    """
    741     Compare gamma_1(B) and gamma_2(B).
     847    Compare gamma_1(S1) and gamma_2(S2).
    742848   
    743     Return return -1 if gamma_1(B) < gamma_2(B), 0 if gamma_1(B) == gamma_2(B),
    744     1 if gamma_1(B) > gamma_2(B).  (Just like the python \code{cmp}) function.
     849    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
     850    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
    745851   
    746852    INPUT:
    747853    gamma_1, gamma_2 -- list permutations (inverse)
    748     S -- a binary code struct object
     854    S1, S2 -- a binary code struct object
    749855   
    750856    """
    751857    cdef int side=0, i, start, end, n_one_1, n_one_2, cur_col
    752858    cdef int where_0, where_1
    753     cdef NonlinearBinaryCodeStruct BCS = <NonlinearBinaryCodeStruct> S
    754     cdef bitset_s *B_1_0 = BCS.scratch_bitsets                   # nwords of len degree
    755     cdef bitset_s *B_1_1 = &BCS.scratch_bitsets[BCS.nwords]      # nwords of len degree
    756     cdef bitset_s *B_2_0 = &BCS.scratch_bitsets[2*BCS.nwords]    # nwords of len degree
    757     cdef bitset_s *B_2_1 = &BCS.scratch_bitsets[3*BCS.nwords]    # nwords of len degree
    758     cdef bitset_s *dividers = &BCS.scratch_bitsets[4*BCS.nwords] # 1 of len nwords
     859    cdef NonlinearBinaryCodeStruct BCS1 = <NonlinearBinaryCodeStruct> S1
     860    cdef NonlinearBinaryCodeStruct BCS2 = <NonlinearBinaryCodeStruct> S2
     861    cdef bitset_s *B_1_0 = BCS1.scratch_bitsets                   # nwords of len degree
     862    cdef bitset_s *B_1_1 = &BCS1.scratch_bitsets[BCS1.nwords]      # nwords of len degree
     863    cdef bitset_s *B_2_0 = &BCS1.scratch_bitsets[2*BCS1.nwords]    # nwords of len degree
     864    cdef bitset_s *B_2_1 = &BCS1.scratch_bitsets[3*BCS1.nwords]    # nwords of len degree
     865    cdef bitset_s *dividers = &BCS1.scratch_bitsets[4*BCS1.nwords] # 1 of len nwords
    759866    cdef bitset_s *B_1_this, *B_1_other, *B_2_this, *B_2_other
    760     for i from 0 <= i < BCS.nwords:
    761         bitset_copy(&B_1_0[i], &BCS.words[i])
    762         bitset_copy(&B_2_0[i], &BCS.words[i])
     867    for i from 0 <= i < BCS1.nwords:
     868        bitset_copy(&B_1_0[i], &BCS1.words[i])
     869        bitset_copy(&B_2_0[i], &BCS2.words[i])
    763870    bitset_zero(dividers)
    764     bitset_set(dividers, BCS.nwords-1)
     871    bitset_set(dividers, BCS1.nwords-1)
    765872   
    766     for cur_col from 0 <= cur_col < BCS.degree:
     873    for cur_col from 0 <= cur_col < BCS1.degree:
    767874        if side == 0:
    768875            B_1_this  = B_1_0
    769876            B_1_other = B_1_1
     
    776883            B_2_other = B_2_0
    777884        side ^= 1
    778885        start = 0
    779         while start < BCS.nwords:
     886        while start < BCS1.nwords:
    780887            end = start
    781888            while not bitset_check(dividers, end):
    782889                end += 1
     
    9601067def random_tests(t=10.0, n_max=50, k_max=6, nwords_max=200, perms_per_code=10, density_range=(.1,.9)):
    9611068    """
    9621069    Tests to make sure that C(gamma(B)) == C(B) for random permutations gamma
    963     and random codes B.
     1070    and random codes B, and that is_isomorphic returns an isomorphism.
    9641071
    9651072    INPUT:
    9661073    t -- run tests for approximately this many seconds
     
    9771084
    9781085    For each code B (B2) generated, we uniformly generate perms_per_code random
    9791086    permutations and verify that the canonical labels of B and the image of B
    980     under the generated permutation are equal.
     1087    under the generated permutation are equal, and check that the double coset
     1088    function returns an isomorphism.
    9811089   
    9821090    DOCTEST:
    9831091        sage: import sage.groups.perm_gps.partn_ref.refinement_binary
     
    9921100    from sage.combinat.permutation import Permutations
    9931101    from sage.matrix.constructor import random_matrix, matrix
    9941102    from sage.rings.finite_field import FiniteField as GF
    995     cdef int h, i, j, n, k, num_tests = 0, num_codes = 0, passed = 1
     1103    cdef int h, i, j, n, k, num_tests = 0, num_codes = 0
    9961104    cdef LinearBinaryCodeStruct B, C
    9971105    cdef NonlinearBinaryCodeStruct B_n, C_n
    9981106    t_0 = walltime()
     
    10401148                    B_n_M[j,B_n_relab[h]] = bitset_check(&B_n.words[j], h)
    10411149                    C_n_M[j,C_n_relab[h]] = bitset_check(&C_n.words[j], h)
    10421150            if B_M.row_space() != C_M.row_space():
    1043                 print "B:"
     1151                print "can_lab error -- B:"
    10441152                for j from 0 <= j < B.dimension:
    10451153                    print bitset_string(&B.basis[j])
    10461154                print perm
    10471155                return
    10481156            if sorted(B_n_M.rows()) != sorted(C_n_M.rows()):
    1049                 print "B_n:"
     1157                print "can_lab error -- B_n:"
    10501158                for j from 0 <= j < B_n.nwords:
    10511159                    print bitset_string(&B_n.words[j])
    10521160                print perm
    10531161                return
     1162            isom = B.is_isomorphic(C)
     1163            if not isom:
     1164                print "isom -- B:"
     1165                for j from 0 <= j < B.dimension:
     1166                    print bitset_string(&B.basis[j])
     1167                print perm
     1168                print isom
     1169                return
     1170            isom = B_n.is_isomorphic(C_n)
     1171            if not isom:
     1172                print "isom -- B_n:"
     1173                for j from 0 <= j < B_n.nwords:
     1174                    print bitset_string(&B_n.words[j])
     1175                print perm
     1176                print isom
     1177                return
    10541178               
    1055         num_tests += perms_per_code
     1179        num_tests += 4*perms_per_code
    10561180        num_codes += 2
    1057     if passed:
    1058         print "All passed: %d random tests on %d codes."%(num_tests, num_codes)
     1181       
     1182    print "All passed: %d random tests on %d codes."%(num_tests, num_codes)
    10591183
    10601184
  • sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_graphs.pxd
    a b  
    88
    99include '../../../ext/cdefs.pxi'
    1010include '../../../ext/stdsage.pxi'
     11include 'data_structures_pxd.pxi' # includes bitsets
    1112
    1213from sage.graphs.base.c_graph cimport CGraph
    1314from sage.graphs.base.sparse_graph cimport SparseGraph
    1415from sage.graphs.base.dense_graph cimport DenseGraph
    1516from sage.rings.integer cimport Integer
    16 from 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
     17from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     18from double_coset cimport double_coset
    1719
    1820cdef class GraphStruct:
    1921    cdef CGraph G
     
    2224    cdef int *scratch # length 3n+1
    2325
    2426cdef int refine_by_degree(PartitionStack *, object, int *, int)
    25 cdef int compare_graphs(int *, int *, object)
     27cdef int compare_graphs(int *, int *, object, object)
    2628cdef bint all_children_are_equivalent(PartitionStack *, object)
    2729cdef inline int degree(PartitionStack *, CGraph, int, int, bint)
    2830cdef inline int sort_by_function(PartitionStack *, int, int *)
  • sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_graphs.pyx
    a b  
    1717# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
    1818#                         http://www.gnu.org/licenses/
    1919#*****************************************************************************
     20
     21include 'data_structures_pyx.pxi' # includes bitsets
     22
     23def isomorphic(G1, G2, partition, ordering2, dig, use_indicator_function, sparse=False):
     24    """
     25    Tests whether two graphs are isomorphic.
     26   
     27    sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import isomorphic
     28
     29    sage: G = Graph(2)
     30    sage: H = Graph(2)
     31    sage: isomorphic(G, H, [[0,1]], [0,1], 0, 1)
     32    [0, 1]
     33    sage: isomorphic(G, H, [[0,1]], [0,1], 0, 1)
     34    [0, 1]
     35    sage: isomorphic(G, H, [[0],[1]], [0,1], 0, 1)
     36    [0, 1]
     37    sage: isomorphic(G, H, [[0],[1]], [1,0], 0, 1)
     38    [1, 0]
     39
     40    sage: G = Graph(3)
     41    sage: H = Graph(3)
     42    sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1)
     43    [0, 1, 2]
     44    sage: G.add_edge(0,1)
     45    sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1)
     46    False
     47    sage: H.add_edge(1,2)
     48    sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1)
     49    [1, 2, 0]
     50
     51    """
     52    cdef int **part
     53    cdef int *output, *ordering
     54    cdef CGraph G
     55    cdef GraphStruct GS1 = GraphStruct()
     56    cdef GraphStruct GS2 = GraphStruct()
     57    cdef GraphStruct GS
     58    cdef int i, j, n = -1
     59
     60    from sage.graphs.graph import GenericGraph, Graph, DiGraph
     61    for G_in in [G1, G2]:
     62        if G_in is G1:
     63            GS = GS1
     64        else:
     65            GS = GS2
     66        if isinstance(G_in, GenericGraph):
     67            if n == -1:
     68                n = G_in.num_verts()
     69            elif n != G_in.num_verts():
     70                return False
     71            G_in = G_in.copy()
     72            if G_in.vertices() != range(n):
     73                to = G_in.relabel(return_map=True)
     74                frm = {}
     75                for v in to.iterkeys():
     76                    frm[to[v]] = v
     77                partition = [[to[v] for v in cell] for cell in partition]
     78            else:
     79                to = range(n)
     80                frm = to
     81            if sparse:
     82                G = SparseGraph(n)
     83            else:
     84                G = DenseGraph(n)
     85            if G_in.is_directed():
     86                for i from 0 <= i < n:
     87                    for _,j,_ in G_in.outgoing_edge_iterator(i):
     88                        G.add_arc(i,j)
     89            else:
     90                for i from 0 <= i < n:
     91                    for _,j,_ in G_in.edge_iterator(i):
     92                        if j <= i:
     93                            G.add_arc(i,j)
     94                            G.add_arc(j,i)
     95        elif isinstance(G_in, CGraph):
     96            G = <CGraph> G_in
     97            if n == -1:
     98                n = G.num_verts
     99            elif n != G.num_verts:
     100                return False
     101            to = range(n)
     102            frm = to
     103        else:
     104            raise TypeError("G must be a Sage graph.")
     105        GS.G = G
     106        GS.directed = 1 if dig else 0
     107        GS.use_indicator = 1 if use_indicator_function else 0
     108
     109    if n == 0:
     110        return {}
     111   
     112    part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     113    ordering = <int *> sage_malloc(n * sizeof(int))
     114    if part is NULL or ordering is NULL:
     115        if part is not NULL: sage_free(part)
     116        if ordering is not NULL: sage_free(ordering)
     117        raise MemoryError
     118    for i from 0 <= i < len(partition):
     119        part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
     120        if part[i] is NULL:
     121            for j from 0 <= j < i:
     122                sage_free(part[j])
     123            sage_free(part)
     124            raise MemoryError
     125        for j from 0 <= j < len(partition[i]):
     126            part[i][j] = partition[i][j]
     127        part[i][len(partition[i])] = -1
     128    part[len(partition)] = NULL
     129    for i from 0 <= i < n:
     130        ordering[i] = ordering2[i]
     131
     132    GS1.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
     133    GS2.scratch = <int *> sage_malloc((3*n+1) * sizeof(int))
     134    if GS1.scratch is NULL or GS2.scratch is NULL:
     135        if GS1.scratch is not NULL: sage_free(GS1.scratch)
     136        if GS2.scratch is not NULL: sage_free(GS2.scratch)
     137        for j from 0 <= j < len(partition):
     138            sage_free(part[j])
     139        sage_free(part)
     140        raise MemoryError
     141
     142    output = double_coset(GS1, GS2, part, ordering, n, &all_children_are_equivalent, &refine_by_degree, &compare_graphs)
     143
     144    for i from 0 <= i < len(partition):
     145        sage_free(part[i])
     146    sage_free(part)
     147    sage_free(ordering)
     148    sage_free(GS1.scratch)
     149    sage_free(GS2.scratch)
     150   
     151    if output is NULL:
     152        return False
     153    else:
     154        output_py = [output[i] for i from 0 <= i < n]
     155        sage_free(output)
     156        # TODO: figure out frm, to stuff to relabel for consistency with input
     157        return output_py
    20158
    21159def search_tree(G_in, partition, lab=True, dig=False, dict_rep=False, certify=False,
    22160                    verbosity=0, use_indicator_function=True, sparse=False,
     
    200338        sage: st(G, [range(G.num_verts())])[1] == st(H, [range(H.num_verts())])[1]
    201339        True
    202340
     341        sage: from sage.graphs.graph import graph_isom_equivalent_non_multi_graph
     342        sage: G = Graph(multiedges=True, implementation='networkx')
     343        sage: G.add_edge(('a', 'b'))
     344        sage: G.add_edge(('a', 'b'))
     345        sage: G.add_edge(('a', 'b'))
     346        sage: G, Pi = graph_isom_equivalent_non_multi_graph(G, [['a','b']])
     347        sage: s,b = st(G, Pi, lab=False, dict_rep=True)
     348        sage: sorted(b.items())
     349        [(('o', 'a'), 5), (('o', 'b'), 1), (('x', 0), 2), (('x', 1), 3), (('x', 2), 4)]
     350
     351        sage: st(Graph(':Dkw'), [range(5)], lab=False, dig=True)
     352        [[4, 1, 2, 3, 0], [0, 2, 1, 3, 4]]
     353
    203354    """
    204355    cdef CGraph G
    205356    cdef int i, j, n
     
    277428    if dict_rep:
    278429        ddd = {}
    279430        for v in frm.iterkeys():
    280             if frm[v] != 0:
    281                 ddd[v] = frm[v]
    282             else:
    283                 ddd[v] = n
     431            ddd[frm[v]] = v if v != 0 else n
    284432        return_tuple.append(ddd)
    285433    if lab:
    286434        if isinstance(G_in, GenericGraph):
     
    348496    cdef int current_cell, i, r
    349497    cdef int first_largest_subcell
    350498    cdef int invariant = 1
     499    cdef int max_degree
    351500    cdef int *degrees = GS.scratch # length 3n+1
    352501    cdef bint necessary_to_split_cell
    353502    cdef int against_index
     
    358507            invariant += 50
    359508            i = current_cell
    360509            necessary_to_split_cell = 0
     510            max_degree = 0
    361511            while 1:
    362512                degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 0)
    363513                if degrees[i-current_cell] != degrees[0]:
    364514                    necessary_to_split_cell = 1
     515                if degrees[i-current_cell] > max_degree:
     516                    max_degree = degrees[i-current_cell]
    365517                i += 1
    366518                if PS.levels[i-1] <= PS.depth:
    367519                    break
     
    369521            if necessary_to_split_cell:
    370522                invariant += 10
    371523                first_largest_subcell = sort_by_function(PS, current_cell, degrees)
    372                 invariant += first_largest_subcell
     524                invariant += first_largest_subcell + max_degree
    373525                against_index = current_cell_against
    374526                while against_index < ctrb_len:
    375527                    if cells_to_refine_by[against_index] == current_cell:
     
    385537                    r += 1
    386538                    if r >= i:
    387539                        break
    388                 invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 0)
    389540                invariant += (i - current_cell)
    390541            current_cell = i
    391542        if GS.directed:
     
    396547                invariant += 20
    397548                i = current_cell
    398549                necessary_to_split_cell = 0
     550                max_degree = 0
    399551                while 1:
    400552                    degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 1)
    401553                    if degrees[i-current_cell] != degrees[0]:
    402554                        necessary_to_split_cell = 1
     555                    if degrees[i-current_cell] > max_degree:
     556                        max_degree = degrees[i-current_cell]
    403557                    i += 1
    404558                    if PS.levels[i-1] <= PS.depth:
    405559                        break
     
    407561                if necessary_to_split_cell:
    408562                    invariant += 7
    409563                    first_largest_subcell = sort_by_function(PS, current_cell, degrees)
    410                     invariant += first_largest_subcell
     564                    invariant += first_largest_subcell + max_degree
    411565                    against_index = current_cell_against
    412566                    while against_index < ctrb_len:
    413567                        if cells_to_refine_by[against_index] == current_cell:
     
    425579                        r += 1
    426580                        if r >= i:
    427581                            break
    428                     invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 1)
    429582                    invariant += (i - current_cell)
    430583                current_cell = i
    431584        current_cell_against += 1
     
    434587    else:
    435588        return 0
    436589   
    437 cdef int compare_graphs(int *gamma_1, int *gamma_2, object S):
     590cdef int compare_graphs(int *gamma_1, int *gamma_2, object S1, object S2):
    438591    r"""
    439     Compare gamma_1(G) and gamma_2(G).
     592    Compare gamma_1(S1) and gamma_2(S2).
    440593
    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
     594    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) ==
     595    gamma_2(S2), 1 if gamma_1(S1) > gamma_2(S2).  (Just like the python
    443596    \code{cmp}) function.
    444597   
    445598    INPUT:
    446599    gamma_1, gamma_2 -- list permutations (inverse)
    447     S -- a graph struct object
     600    S1, S2 -- graph struct objects
    448601   
    449602    """
    450603    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]):
     604    cdef GraphStruct GS1 = <GraphStruct> S1
     605    cdef GraphStruct GS2 = <GraphStruct> S2
     606    cdef CGraph G1 = GS1.G
     607    cdef CGraph G2 = GS2.G
     608    for i from 0 <= i < G1.num_verts:
     609        for j from 0 <= j < G1.num_verts:
     610            if G1.has_arc_unsafe(gamma_1[i], gamma_1[j]):
     611                if not G2.has_arc_unsafe(gamma_2[i], gamma_2[j]):
    457612                    return 1
    458             elif G.has_arc_unsafe(gamma_2[i], gamma_2[j]):
     613            elif G2.has_arc_unsafe(gamma_2[i], gamma_2[j]):
    459614                return -1
    460615    return 0
    461616
     
    472627    S -- a graph struct object
    473628    """
    474629    cdef GraphStruct GS = <GraphStruct> S
     630    if GS.directed:
     631        return 0
    475632    cdef CGraph G = GS.G
    476633    cdef int i, n = PS.degree
    477634    cdef bint in_cell = 0
     
    622779def random_tests(t=10.0, n_max=60, perms_per_graph=10):
    623780    """
    624781    Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma
    625     and random graphs G.
     782    and random graphs G, and that isomorphic returns an isomorphism.
    626783
    627784    INPUT:
    628785    t -- run tests for approximately this many seconds
     
    637794
    638795    For each graph G generated, we uniformly generate perms_per_graph random
    639796    permutations and verify that the canonical labels of G and the image of G
    640     under the generated permutation are equal.
     797    under the generated permutation are equal, and that the isomorphic function
     798    returns an isomorphism.
    641799       
    642800    TESTS:
    643801        sage: import sage.groups.perm_gps.partn_ref.refinement_graphs
     
    651809    from sage.misc.prandom import random, randint
    652810    from sage.graphs.graph_generators import GraphGenerators, DiGraphGenerators
    653811    from sage.combinat.permutation import Permutations
    654     cdef int i, j, num_tests = 0, num_graphs = 0, passed = 1
     812    cdef int i, j, num_tests = 0, num_graphs = 0
    655813    GG = GraphGenerators()
    656814    DGG = DiGraphGenerators()
    657815    t_0 = walltime()
     
    663821        G = GG.RandomGNP(n, p)
    664822        H = G.copy()
    665823        for i from 0 <= i < perms_per_graph:
     824            G = H.copy()
    666825            G1 = search_tree(G, [G.vertices()])[1]
    667826            perm = list(S.random_element())
    668827            perm = [perm[j]-1 for j from 0 <= j < n]
    669828            G.relabel(perm)
    670829            G2 = search_tree(G, [G.vertices()])[1]
    671830            if G1 != G2:
    672                 print "FAILURE: graph6-"
     831                print "search_tree FAILURE: graph6-"
    673832                print H.graph6_string()
    674833                print perm
    675                 passed = 0
     834                return
     835            isom = isomorphic(G, H, [range(n)], range(n), 0, 1)
     836            if not isom or G.relabel(isom, inplace=False) != H:
     837                print "isom FAILURE: graph6-"
     838                print H.graph6_string()
     839                print perm
     840                return
    676841
    677842        D = DGG.RandomDirectedGNP(n, p)
    678843        D.loops(True)
     
    681846                D.add_edge(i,i)
    682847        E = D.copy()
    683848        for i from 0 <= i < perms_per_graph:
     849            D = E.copy()
    684850            D1 = search_tree(D, [D.vertices()], dig=True)[1]
    685851            perm = list(S.random_element())
    686852            perm = [perm[j]-1 for j from 0 <= j < n]
    687853            D.relabel(perm)
    688854            D2 = search_tree(D, [D.vertices()], dig=True)[1]
    689855            if D1 != D2:
    690                 print "FAILURE: dig6-"
     856                print "search_tree FAILURE: dig6-"
    691857                print E.dig6_string()
    692858                print perm
    693                 passed = 0
    694         num_tests += 2*perms_per_graph
     859                return
     860            isom = isomorphic(D, E, [range(n)], range(n), 1, 1)
     861            if not isom or D.relabel(isom, inplace=False) != E:
     862                print "isom FAILURE: dig6-"
     863                print E.dig6_string()
     864                print perm
     865                print isom
     866                return
     867        num_tests += 4*perms_per_graph
    695868        num_graphs += 2
    696     if passed:
    697         print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs)
     869    print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs)
    698870
    699871
  • sage/groups/perm_gps/partn_ref/refinement_matrices.pxd

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_matrices.pxd
    a b  
    88
    99include '../../../ext/cdefs.pxi'
    1010include '../../../ext/stdsage.pxi'
    11 include '../../../misc/bitset_pxd.pxi'
     11include 'data_structures_pxd.pxi' # includes bitsets
    1212
    1313from sage.rings.integer cimport Integer
    14 from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport PartitionStack, PS_new, PS_dealloc, PS_copy_from_to, get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
    15 from sage.groups.perm_gps.partn_ref.refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree
    16 from sage.groups.perm_gps.partn_ref.refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent
     14from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return
     15from refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree
     16from refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent
     17from double_coset cimport double_coset
    1718
    1819cdef class MatrixStruct:
    1920    cdef list symbol_structs
     
    2627    cdef aut_gp_and_can_lab_return *output
    2728
    2829cdef int refine_matrix(PartitionStack *, object, int *, int)
    29 cdef int compare_matrices(int *, int *, object)
     30cdef int compare_matrices(int *, int *, object, object)
    3031cdef bint all_matrix_children_are_equivalent(PartitionStack *, object)
    3132
  • sage/groups/perm_gps/partn_ref/refinement_matrices.pyx

    diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_matrices.pyx
    a b  
    2222#                         http://www.gnu.org/licenses/
    2323#*****************************************************************************
    2424
    25 include '../../../misc/bitset.pxi'
     25include 'data_structures_pyx.pxi' # includes bitsets
     26
    2627from sage.misc.misc import uniq
    2728from sage.matrix.constructor import Matrix
    2829
     
    127128
    128129        """
    129130        cdef int **part, i, j
     131        cdef NonlinearBinaryCodeStruct S_temp
     132        for i from 0 <= i < self.nsymbols:
     133            S_temp = <NonlinearBinaryCodeStruct> self.symbol_structs[i]
     134            S_temp.first_time = 1
     135
    130136        if partition is None:
    131137            partition = [range(self.degree)]
    132138        part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     
    195201            self.run()
    196202        return [self.output.relabeling[i] for i from 0 <= i < self.degree]
    197203
     204    def is_isomorphic(self, MatrixStruct other):
     205        """
     206        Calculate whether self is isomorphic to other.
     207       
     208        EXAMPLES:
     209            sage: from sage.groups.perm_gps.partn_ref.refinement_matrices import MatrixStruct
     210            sage: M = MatrixStruct(Matrix(GF(11), [[1,2,3,0,0,0],[0,0,0,1,2,3]]))
     211            sage: N = MatrixStruct(Matrix(GF(11), [[0,1,0,2,0,3],[1,0,2,0,3,0]]))
     212            sage: M.is_isomorphic(N)
     213            [0, 2, 4, 1, 3, 5]
     214           
     215        """
     216       
     217        cdef int **part, i, j
     218        cdef int *output, *ordering
     219        cdef NonlinearBinaryCodeStruct S_temp
     220        for i from 0 <= i < self.nsymbols:
     221            S_temp = self.symbol_structs[i]
     222            S_temp.first_time = 1
     223            S_temp = other.symbol_structs[i]
     224            S_temp.first_time = 1
     225        partition = [range(self.degree)]
     226        part = <int **> sage_malloc((len(partition)+1) * sizeof(int *))
     227        ordering = <int *> sage_malloc(self.degree * sizeof(int))
     228        if part is NULL or ordering is NULL:
     229            if part is not NULL: sage_free(part)
     230            if ordering is not NULL: sage_free(ordering)
     231            raise MemoryError
     232        for i from 0 <= i < len(partition):
     233            part[i] = <int *> sage_malloc((len(partition[i])+1) * sizeof(int))
     234            if part[i] is NULL:
     235                for j from 0 <= j < i:
     236                    sage_free(part[j])
     237                sage_free(part)
     238                raise MemoryError
     239            for j from 0 <= j < len(partition[i]):
     240                part[i][j] = partition[i][j]
     241            part[i][len(partition[i])] = -1
     242        part[len(partition)] = NULL
     243        for i from 0 <= i < self.degree:
     244            ordering[i] = i
     245       
     246        output = double_coset(self, other, part, ordering, self.degree, &all_matrix_children_are_equivalent, &refine_matrix, &compare_matrices)
     247       
     248        for i from 0 <= i < len(partition):
     249            sage_free(part[i])
     250        sage_free(part)
     251        sage_free(ordering)
     252
     253        if output is NULL:
     254            return False
     255        else:
     256            output_py = [output[i] for i from 0 <= i < self.degree]
     257            sage_free(output)
     258            return output_py
     259   
    198260cdef int refine_matrix(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len):
    199261    cdef MatrixStruct M = <MatrixStruct> S
    200262    cdef int i, temp_inv, invariant = 1
     
    209271            changed = 0
    210272    return invariant
    211273   
    212 cdef int compare_matrices(int *gamma_1, int *gamma_2, object S):
    213     cdef MatrixStruct MS = <MatrixStruct> S
    214     M = MS.matrix
     274cdef int compare_matrices(int *gamma_1, int *gamma_2, object S1, object S2):
     275    cdef MatrixStruct MS1 = <MatrixStruct> S1
     276    cdef MatrixStruct MS2 = <MatrixStruct> S2
     277    M1 = MS1.matrix
     278    M2 = MS2.matrix
    215279    cdef int i
    216     M1 = Matrix(M.base_ring(), M.nrows(), M.ncols(), sparse=M.is_sparse())
    217     M2 = Matrix(M.base_ring(), M.nrows(), M.ncols(), sparse=M.is_sparse())
    218     for i from 0 <= i < M.ncols():
    219         M1.set_column(i, M.column(gamma_1[i]))
    220         M2.set_column(i, M.column(gamma_2[i]))
    221     return cmp(sorted(M1.rows()), sorted(M2.rows()))
     280    MM1 = Matrix(M1.base_ring(), M1.nrows(), M1.ncols(), sparse=M1.is_sparse())
     281    MM2 = Matrix(M2.base_ring(), M2.nrows(), M2.ncols(), sparse=M2.is_sparse())
     282    for i from 0 <= i < M1.ncols():
     283        MM1.set_column(i, M1.column(gamma_1[i]))
     284        MM2.set_column(i, M2.column(gamma_2[i]))
     285    return cmp(sorted(MM1.rows()), sorted(MM2.rows()))
    222286
    223287cdef bint all_matrix_children_are_equivalent(PartitionStack *PS, object S):
    224288    cdef MatrixStruct M = <MatrixStruct> S
     
    230294def random_tests(t=10.0, nrows_max=50, ncols_max=50, nsymbols_max=20, perms_per_matrix=10, density_range=(.1,.9)):
    231295    """
    232296    Tests to make sure that C(gamma(M)) == C(M) for random permutations gamma
    233     and random matrices M.
     297    and random matrices M, and that M.is_isomorphic(gamma(M)) returns an
     298    isomorphism.
    234299
    235300    INPUT:
    236301    t -- run tests for approximately this many seconds
     
    247312
    248313    For each matrix M generated, we uniformly generate perms_per_matrix random
    249314    permutations and verify that the canonical labels of M and the image of M
    250     under the generated permutation are equal.
     315    under the generated permutation are equal, and that the isomorphism is
     316    discovered by the double coset function.
    251317   
    252318    DOCTEST:
    253319        sage: import sage.groups.perm_gps.partn_ref.refinement_matrices
     
    263329    from sage.matrix.constructor import random_matrix, matrix
    264330    from sage.rings.finite_field import FiniteField as GF
    265331    from sage.rings.arith import next_prime
    266     cdef int h, i, j, nrows, k, num_tests = 0, num_matrices = 0, passed = 1
     332    cdef int h, i, j, nrows, k, num_tests = 0, num_matrices = 0
    267333    cdef MatrixStruct M, N
    268334    t_0 = walltime()
    269335    while walltime(t_0) < t:
     
    300366
    301367            if M_C != N_C:
    302368                print "M:"
    303                 print M
     369                print M.matrix.str()
    304370                print "perm:"
    305371                print perm
     372                return
     373
     374            isom = M.is_isomorphic(N)
     375            if not isom:
     376                print "isom FAILURE: M:"
     377                print M.matrix.str()
     378                print "isom FAILURE: N:"
     379                print N.matrix.str()
    306380                return
    307381               
    308382        num_tests += perms_per_matrix
    309383        num_matrices += 2
    310     if passed:
    311         print "All passed: %d random tests on %d matrices."%(num_tests, num_matrices)
     384    print "All passed: %d random tests on %d matrices."%(num_tests, num_matrices)
    312385
    313386
    314387
  • setup.py

    diff -r 0d158c58d1c2 -r 2df5cad14396 setup.py
    a b  
    717717
    718718    Extension('sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label',
    719719              sources = ['sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx']), \
     720
     721    Extension('sage.groups.perm_gps.partn_ref.double_coset',
     722              sources = ['sage/groups/perm_gps/partn_ref/double_coset.pyx']), \
    720723
    721724    Extension('sage.groups.perm_gps.partn_ref.refinement_graphs',
    722725              sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx']), \