# HG changeset patch # User Robert L. Miller # 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/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd Mon Sep 15 17:34:52 2008 -0700 @@ -8,55 +8,9 @@ include '../../../ext/cdefs.pxi' include '../../../ext/stdsage.pxi' -include '../../../misc/bitset_pxd.pxi' +include 'data_structures_pxd.pxi' # includes bitsets from sage.rings.integer cimport Integer - -cdef struct OrbitPartition: - # Disjoint set data structure for representing the orbits of the generators - # so far found. Also keeps track of the minimum elements of the cells and - # their sizes. - int degree - int *parent - int *rank - int *mcr # minimum cell representatives - only valid at the root of a cell - int *size # also only valid at the root of a cell - -cdef inline OrbitPartition *OP_new(int) -cdef inline int OP_dealloc(OrbitPartition *) -cdef inline int OP_find(OrbitPartition *, int) -cdef inline int OP_join(OrbitPartition *, int, int) -cdef inline int OP_merge_list_perm(OrbitPartition *, int *) - -cdef struct PartitionStack: - # Representation of a node of the search tree. A sequence of partitions of - # length depth + 1, each of which is finer than the last. Partition k is - # represented as PS.entries in order, broken up immediately after each - # entry of levels which is at most k. - int *entries - int *levels - int depth - int degree - -cdef inline PartitionStack *PS_new(int, bint) -cdef inline PartitionStack *PS_copy(PartitionStack *) -cdef inline int PS_copy_from_to(PartitionStack *, PartitionStack *) -cdef inline PartitionStack *PS_from_list(object) -cdef inline int PS_dealloc(PartitionStack *) -cdef inline int PS_print(PartitionStack *) -cdef inline int PS_print_partition(PartitionStack *, int) -cdef inline bint PS_is_discrete(PartitionStack *) -cdef inline int PS_num_cells(PartitionStack *) -cdef inline int PS_move_min_to_front(PartitionStack *, int, int) -cdef inline bint PS_is_mcr(PartitionStack *, int) -cdef inline bint PS_is_fixed(PartitionStack *, int) -cdef inline int PS_clear(PartitionStack *) -cdef inline int PS_split_point(PartitionStack *, int) -cdef inline int PS_first_smallest(PartitionStack *, bitset_t) -cdef inline int PS_get_perm_from(PartitionStack *, PartitionStack *, int *) - -cdef inline int split_point_and_refine(PartitionStack *, int, object, - int (*)(PartitionStack *, object, int *, int), int *) cdef struct aut_gp_and_can_lab_return: int *generators @@ -69,7 +23,7 @@ cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab( object, int **, int, bint (*)(PartitionStack *, object), int (*)(PartitionStack *, object, int *, int), - int (*)(int *, int *, object), bint, bint, bint) + int (*)(int *, int *, object, object), bint, bint, bint) diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx --- a/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx Mon Sep 15 17:34:52 2008 -0700 @@ -58,11 +58,11 @@ Signature: - \code{int compare_structures(int *gamma_1, int *gamma_2, object S)} + \code{int compare_structures(int *gamma_1, int *gamma_2, object S1, object S2)} This function must implement a total ordering on the set of objects of fixed - order. Return -1 if \code{gamma_1(S) < gamma_2(S)}, 0 if - \code{gamma_1(S) == gamma_2(S)}, 1 if \code{gamma_1(S) > gamma_2(S)}. + order. Return -1 if \code{gamma_1(S1) < gamma_2(S2)}, 0 if + \code{gamma_1(S1) == gamma_2(S2)}, 1 if \code{gamma_1(S1) > gamma_2(S2)}. C. \code{all_children_are_equivalent}: @@ -93,600 +93,14 @@ # http://www.gnu.org/licenses/ #***************************************************************************** -include '../../../misc/bitset.pxi' - -cdef inline int min(int a, int b): - if a < b: return a - else: return b - -cdef inline int max(int a, int b): - if a > b: return a - else: return b +include 'data_structures_pyx.pxi' # includes bitsets cdef inline int cmp(int a, int b): if a < b: return -1 elif a == b: return 0 else: return 1 -# OrbitPartitions - -cdef inline OrbitPartition *OP_new(int n): - """ - Allocate and return a pointer to a new OrbitPartition of degree n. Returns a - null pointer in the case of an allocation failure. - """ - cdef int i - cdef OrbitPartition *OP = \ - sage_malloc(sizeof(OrbitPartition)) - if OP is NULL: - return OP - OP.degree = n - OP.parent = sage_malloc( n * sizeof(int) ) - OP.rank = sage_malloc( n * sizeof(int) ) - OP.mcr = sage_malloc( n * sizeof(int) ) - OP.size = sage_malloc( n * sizeof(int) ) - if OP.parent is NULL or OP.rank is NULL or \ - OP.mcr is NULL or OP.size is NULL: - if OP.parent is not NULL: sage_free(OP.parent) - if OP.rank is not NULL: sage_free(OP.rank) - if OP.mcr is not NULL: sage_free(OP.mcr) - if OP.size is not NULL: sage_free(OP.size) - return NULL - for i from 0 <= i < n: - OP.parent[i] = i - OP.rank[i] = 0 - OP.mcr[i] = i - OP.size[i] = 1 - return OP - -cdef inline int OP_dealloc(OrbitPartition *OP): - sage_free(OP.parent) - sage_free(OP.rank) - sage_free(OP.mcr) - sage_free(OP.size) - sage_free(OP) - return 0 - -cdef inline int OP_find(OrbitPartition *OP, int n): - """ - Report the representative ("root") of the cell which contains n. - """ - if OP.parent[n] == n: - return n - else: - OP.parent[n] = OP_find(OP, OP.parent[n]) - return OP.parent[n] - -cdef inline int OP_join(OrbitPartition *OP, int m, int n): - """ - Join the cells containing m and n, if they are different. - """ - cdef int m_root = OP_find(OP, m) - cdef int n_root = OP_find(OP, n) - if OP.rank[m_root] > OP.rank[n_root]: - OP.parent[n_root] = m_root - OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root]) - OP.size[m_root] += OP.size[n_root] - elif OP.rank[m_root] < OP.rank[n_root]: - OP.parent[m_root] = n_root - OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root]) - OP.size[n_root] += OP.size[m_root] - elif m_root != n_root: - OP.parent[n_root] = m_root - OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root]) - OP.size[m_root] += OP.size[n_root] - OP.rank[m_root] += 1 - -cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma): - """ - Joins the cells of OP which intersect the same orbit of gamma. - - INPUT: - gamma - an integer array representing i -> gamma[i]. - - OUTPUT: - 1 - something changed - 0 - orbits of gamma all contained in cells of OP - """ - cdef int i, i_root, gamma_i_root, changed = 0 - for i from 0 <= i < OP.degree: - if gamma[i] == i: continue - i_root = OP_find(OP, i) - gamma_i_root = OP_find(OP, gamma[i]) - if i_root != gamma_i_root: - changed = 1 - OP_join(OP, i_root, gamma_i_root) - return changed - -def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]): - """ - Demonstration and testing. - - DOCTEST: - sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \ - ... import OP_represent - sage: OP_represent() - Allocating OrbitPartition... - Allocation passed. - Checking that each element reports itself as its root. - Each element reports itself as its root. - Merging: - Merged 0 and 1. - Merged 2 and 3. - Merged 3 and 4. - Done merging. - Finding: - 0 -> 0, root: size=2, mcr=0, rank=1 - 1 -> 0 - 2 -> 2, root: size=3, mcr=2, rank=1 - 3 -> 2 - 4 -> 2 - 5 -> 5, root: size=1, mcr=5, rank=0 - 6 -> 6, root: size=1, mcr=6, rank=0 - 7 -> 7, root: size=1, mcr=7, rank=0 - 8 -> 8, root: size=1, mcr=8, rank=0 - Allocating array to test merge_perm. - Allocation passed. - Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8] - Done merging. - Finding: - 0 -> 0, root: size=5, mcr=0, rank=2 - 1 -> 0 - 2 -> 0 - 3 -> 0 - 4 -> 0 - 5 -> 5, root: size=3, mcr=5, rank=1 - 6 -> 5 - 7 -> 5 - 8 -> 8, root: size=1, mcr=8, rank=0 - Deallocating OrbitPartition. - Done. - - """ - cdef int i - print "Allocating OrbitPartition..." - cdef OrbitPartition *OP = OP_new(n) - if OP is NULL: - print "Allocation failed!" - return - print "Allocation passed." - print "Checking that each element reports itself as its root." - good = True - for i from 0 <= i < n: - if not OP_find(OP, i) == i: - print "Failed at i = %d!"%i - good = False - if good: print "Each element reports itself as its root." - print "Merging:" - for i,j in merges: - OP_join(OP, i, j) - print "Merged %d and %d."%(i,j) - print "Done merging." - print "Finding:" - for i from 0 <= i < n: - j = OP_find(OP, i) - s = "%d -> %d"%(i, j) - if i == j: - s += ", root: size=%d, mcr=%d, rank=%d"%\ - (OP.size[i], OP.mcr[i], OP.rank[i]) - print s - print "Allocating array to test merge_perm." - cdef int *gamma = sage_malloc( n * sizeof(int) ) - if gamma is NULL: - print "Allocation failed!" - OP_dealloc(OP) - return - print "Allocation passed." - for i from 0 <= i < n: - gamma[i] = perm[i] - print "Merging permutation: %s"%perm - OP_merge_list_perm(OP, gamma) - print "Done merging." - print "Finding:" - for i from 0 <= i < n: - j = OP_find(OP, i) - s = "%d -> %d"%(i, j) - if i == j: - s += ", root: size=%d, mcr=%d, rank=%d"%\ - (OP.size[i], OP.mcr[i], OP.rank[i]) - print s - print "Deallocating OrbitPartition." - sage_free(gamma) - OP_dealloc(OP) - print "Done." - -# PartitionStacks - -cdef inline PartitionStack *PS_new(int n, bint unit_partition): - """ - Allocate and return a pointer to a new PartitionStack of degree n. Returns a - null pointer in the case of an allocation failure. - """ - cdef int i - cdef PartitionStack *PS = \ - sage_malloc(sizeof(PartitionStack)) - if PS is NULL: - return PS - PS.entries = sage_malloc( n * sizeof(int) ) - PS.levels = sage_malloc( n * sizeof(int) ) - if PS.entries is NULL or PS.levels is NULL: - if PS.entries is not NULL: sage_free(PS.entries) - if PS.levels is not NULL: sage_free(PS.levels) - return NULL - PS.depth = 0 - PS.degree = n - if unit_partition: - for i from 0 <= i < n-1: - PS.entries[i] = i - PS.levels[i] = n - PS.entries[n-1] = n-1 - PS.levels[n-1] = -1 - return PS - -cdef inline PartitionStack *PS_copy(PartitionStack *PS): - """ - Allocate and return a pointer to a copy of PartitionStack PS. Returns a null - pointer in the case of an allocation failure. - """ - cdef int i - cdef PartitionStack *PS2 = \ - sage_malloc(sizeof(PartitionStack)) - if PS2 is NULL: - return PS2 - PS2.entries = sage_malloc( PS.degree * sizeof(int) ) - PS2.levels = sage_malloc( PS.degree * sizeof(int) ) - if PS2.entries is NULL or PS2.levels is NULL: - if PS2.entries is not NULL: sage_free(PS2.entries) - if PS2.levels is not NULL: sage_free(PS2.levels) - return NULL - PS_copy_from_to(PS, PS2) - return PS2 - -cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2): - """ - Copy all data from PS to PS2. - """ - PS2.depth = PS.depth - PS2.degree = PS.degree - for i from 0 <= i < PS.degree: - PS2.entries[i] = PS.entries[i] - PS2.levels[i] = PS.levels[i] - -cdef inline PartitionStack *PS_from_list(object L): - """ - Allocate and return a pointer to a PartitionStack representing L. Returns a - null pointer in the case of an allocation failure. - """ - cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0 - for cell from 0 <= cell < num_cells: - n += len(L[cell]) - cdef PartitionStack *PS = PS_new(n, 0) - cell = 0 - if PS is NULL: - return PS - while 1: - cur_len = len(L[cell]) - for i from 0 <= i < cur_len: - PS.entries[cur_start + i] = L[cell][i] - PS.levels[cur_start + i] = n - PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1) - cur_start += cur_len - cell += 1 - if cell == num_cells: - PS.levels[cur_start-1] = -1 - break - PS.levels[cur_start-1] = 0 - return PS - -cdef inline int PS_dealloc(PartitionStack *PS): - sage_free(PS.entries) - sage_free(PS.levels) - sage_free(PS) - return 0 - -cdef inline int PS_print(PartitionStack *PS): - """ - Print a visual representation of PS. - """ - cdef int i - for i from 0 <= i < PS.depth + 1: - PS_print_partition(PS, i) - print - -cdef inline int PS_print_partition(PartitionStack *PS, int k): - """ - Print the partition at depth k. - """ - s = '(' - for i from 0 <= i < PS.degree: - s += str(PS.entries[i]) - if PS.levels[i] <= k: - s += '|' - else: - s += ',' - s = s[:-1] + ')' - print s - -cdef inline bint PS_is_discrete(PartitionStack *PS): - """ - Returns whether the deepest partition consists only of singleton cells. - """ - cdef int i - for i from 0 <= i < PS.degree: - if PS.levels[i] > PS.depth: - return 0 - return 1 - -cdef inline int PS_num_cells(PartitionStack *PS): - """ - Returns the number of cells. - """ - cdef int i, ncells = 0 - for i from 0 <= i < PS.degree: - if PS.levels[i] <= PS.depth: - ncells += 1 - return ncells - -cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end): - """ - Makes sure that the first element of the segment of entries i with - start <= i <= end is minimal. - """ - cdef int i, min_loc = start, minimum = PS.entries[start] - for i from start < i <= end: - if PS.entries[i] < minimum: - min_loc = i - minimum = PS.entries[i] - if min_loc != start: - PS.entries[min_loc] = PS.entries[start] - PS.entries[start] = minimum - -cdef inline bint PS_is_mcr(PartitionStack *PS, int m): - """ - Returns whether PS.elements[m] (not m!) is the smallest element of its cell. - """ - return m == 0 or PS.levels[m-1] <= PS.depth - -cdef inline bint PS_is_fixed(PartitionStack *PS, int m): - """ - Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming - PS_is_mcr(PS, m) is already true. - """ - return PS.levels[m] <= PS.depth - -cdef inline int PS_clear(PartitionStack *PS): - """ - Sets the current partition to the first shallower one, i.e. forgets about - boundaries between cells that are new to the current level. - """ - cdef int i, cur_start = 0 - for i from 0 <= i < PS.degree: - if PS.levels[i] >= PS.depth: - PS.levels[i] += 1 - if PS.levels[i] < PS.depth: - PS_move_min_to_front(PS, cur_start, i) - cur_start = i+1 - -cdef inline int PS_split_point(PartitionStack *PS, int v): - """ - Detaches the point v from the cell it is in, putting the singleton cell - of just v in front. Returns the position where v is now located. - """ - cdef int i = 0, index_of_v - while PS.entries[i] != v: - i += 1 - index_of_v = i - while PS.levels[i] > PS.depth: - i += 1 - if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \ - and PS.levels[index_of_v] > PS.depth: - # if v is first (make sure v is not already alone) - PS_move_min_to_front(PS, index_of_v+1, i) - PS.levels[index_of_v] = PS.depth - return index_of_v - else: - # If v is not at front, i.e. v is not minimal in its cell, - # then move_min_to_front is not necessary since v will swap - # with the first before making its own cell, leaving it at - # the front of the other. - i = index_of_v - while i != 0 and PS.levels[i-1] > PS.depth: - i -= 1 - PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v - PS.entries[i+1] = PS.entries[i] # move the first (min) to second - PS.entries[i] = v # place v first - PS.levels[i] = PS.depth - return i - -cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b): - """ - Find the first occurrence of the smallest cell of size greater than one, - store its entries to b, and return its minimum element. - """ - cdef int i = 0, j = 0, location = 0, n = PS.degree - bitset_zero(b) - while 1: - if PS.levels[i] <= PS.depth: - if i != j and n > i - j + 1: - n = i - j + 1 - location = j - j = i + 1 - if PS.levels[i] == -1: break - i += 1 - # location now points to the beginning of the first, smallest, - # nontrivial cell - i = location - while 1: - bitset_flip(b, PS.entries[i]) - if PS.levels[i] <= PS.depth: break - i += 1 - return PS.entries[location] - -cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma): - """ - Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i] - denotes the entry of the ith cell of the discrete partition PS. - """ - cdef int i = 0 - for i from 0 <= i < PS1.degree: - gamma[PS2.entries[i]] = PS1.entries[i] - -def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]): - """ - Demonstration and testing. - - DOCTEST: - sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \ - ... import PS_represent - sage: PS_represent() - Allocating PartitionStack... - Allocation passed: - (0,1,2,3,4,5,6,7,8,9) - Checking that entries are in order and correct level. - Everything seems in order, deallocating. - Deallocated. - Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]]. - PartitionStack's data: - entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2] - levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1] - depth = 0, degree = 10 - (6|3,4,8,7|1,9,5|0,2) - Checking PS_is_discrete: - False - Checking PS_num_cells: - 4 - Checking PS_is_mcr, min cell reps are: - [6, 3, 1, 0] - Checking PS_is_fixed, fixed elements are: - [6] - Copying PartitionStack: - (6|3,4,8,7|1,9,5|0,2) - Checking for consistency. - Everything is consistent. - Clearing copy: - (0,3,4,8,7,1,9,5,6,2) - Splitting point 6 from original: - 0 - (6|3,4,8,7|1,9,5|0,2) - Splitting point 1 from original: - 5 - (6|3,4,8,7|1|5,9|0,2) - Splitting point 8 from original: - 1 - (6|8|3,4,7|1|5,9|0,2) - Splitting point 2 from original: - 8 - (6|8|3,4,7|1|5,9|2|0) - Getting permutation from PS2->PS: - [6, 1, 0, 8, 3, 9, 2, 7, 4, 5] - Finding first smallest: - Minimal element is 5, bitset is: - 0000010001 - Deallocating PartitionStacks. - Done. - - """ - cdef int i, n = sum([len(cell) for cell in partition]), *gamma - cdef bitset_t b - print "Allocating PartitionStack..." - cdef PartitionStack *PS = PS_new(n, 1), *PS2 - if PS is NULL: - print "Allocation failed!" - return - print "Allocation passed:" - PS_print(PS) - print "Checking that entries are in order and correct level." - good = True - for i from 0 <= i < n-1: - if not (PS.entries[i] == i and PS.levels[i] == n): - print "Failed at i = %d!"%i - print PS.entries[i], PS.levels[i], i, n - good = False - if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1): - print "Failed at i = %d!"%(n-1) - good = False - if not PS.degree == n or not PS.depth == 0: - print "Incorrect degree or depth!" - good = False - if good: print "Everything seems in order, deallocating." - PS_dealloc(PS) - print "Deallocated." - print "Creating PartitionStack from partition %s."%partition - PS = PS_from_list(partition) - print "PartitionStack's data:" - print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n] - print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n] - print "depth = %d, degree = %d"%(PS.depth,PS.degree) - PS_print(PS) - print "Checking PS_is_discrete:" - print "True" if PS_is_discrete(PS) else "False" - print "Checking PS_num_cells:" - print PS_num_cells(PS) - print "Checking PS_is_mcr, min cell reps are:" - L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)] - print L - print "Checking PS_is_fixed, fixed elements are:" - print [PS.entries[l] for l in L if PS_is_fixed(PS, l)] - print "Copying PartitionStack:" - PS2 = PS_copy(PS) - PS_print(PS2) - print "Checking for consistency." - good = True - for i from 0 <= i < n: - if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]: - print "Failed at i = %d!"%i - good = False - if PS.degree != PS2.degree or PS.depth != PS2.depth: - print "Failure with degree or depth!" - good = False - if good: - print "Everything is consistent." - print "Clearing copy:" - PS_clear(PS2) - PS_print(PS2) - for s in splits: - print "Splitting point %d from original:"%s - print PS_split_point(PS, s) - PS_print(PS) - print "Getting permutation from PS2->PS:" - gamma = sage_malloc(n * sizeof(int)) - PS_get_perm_from(PS, PS2, gamma) - print [gamma[i] for i from 0 <= i < n] - sage_free(gamma) - print "Finding first smallest:" - bitset_init(b, n) - i = PS_first_smallest(PS, b) - print "Minimal element is %d, bitset is:"%i - print bitset_string(b) - bitset_clear(b) - print "Deallocating PartitionStacks." - PS_dealloc(PS) - PS_dealloc(PS2) - print "Done." - # Functions - -cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S, - int (*refine_and_return_invariant)\ - (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len), - int *cells_to_refine_by): - """ - Make the partition stack one longer by copying the last partition in the - stack, split off a given point, and refine. Return the invariant given by - the refinement function. - - INPUT: - PS -- the partition stack to refine - v -- the point to split - S -- the structure - refine_and_return_invariant -- the refinement function provided - cells_to_refine_by -- an array, contents ignored - - """ - PS.depth += 1 - PS_clear(PS) - cells_to_refine_by[0] = PS_split_point(PS, v) - return refine_and_return_invariant(PS, S, cells_to_refine_by, 1) cdef bint all_children_are_equivalent_trivial(PartitionStack *PS, object S): return 0 @@ -694,7 +108,7 @@ cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len): return 0 -cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S): +cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S1, object S2): return 0 def test_get_aut_gp_and_can_lab_trivially(int n=6, partition=[[0,1,2],[3,4],[5]], @@ -747,7 +161,7 @@ bint (*all_children_are_equivalent)(PartitionStack *PS, object S), int (*refine_and_return_invariant)\ (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len), - int (*compare_structures)(int *gamma_1, int *gamma_2, object S), + int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2), bint canonical_label, bint base, bint order): """ Traverse the search space for subgroup/canonical label calculation. @@ -779,10 +193,10 @@ int -- returns an invariant under application of arbitrary permutations compare_structures -- pointer to a function INPUT: - gamma_1, gamma_2 -- (list) permutations of the points of S - S -- pointer to the structure + gamma_1, gamma_2 -- (list) permutations of the points of S1 and S2 + S1, S2 -- pointers to the structures OUTPUT: - int -- 0 if gamma_1(S) = gamma_2(S), otherwise -1 or 1 (see docs for cmp), + int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp), such that the set of all structures is well-ordered OUTPUT: pointer to a aut_gp_and_can_lab_return struct @@ -1020,7 +434,8 @@ # Ignore the invariant this time, since we are # creating the root of the search tree. refine_and_return_invariant(current_ps, S, cells_to_refine_by, j) - + PS_move_all_mins_to_front(current_ps) + # Refine down to a discrete partition while not PS_is_discrete(current_ps): if not all_children_are_equivalent(current_ps, S): @@ -1029,6 +444,7 @@ bitset_unset(vertices_have_been_reduced, current_ps.depth) i = current_ps.depth current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by) + PS_move_all_mins_to_front(current_ps) first_indicators[i] = current_indicators[i] if canonical_label: label_indicators[i] = current_indicators[i] @@ -1074,9 +490,8 @@ bitset_flip(vertices_have_been_reduced, current_ps.depth) # Look for a new point to split. i = vertices_determining_current_stack[current_ps.depth] + 1 - while i < n and not bitset_check(vertices_to_split[current_ps.depth], i): - i += 1 - if i < n: + i = bitset_next(vertices_to_split[current_ps.depth], i) + if i != -1: # There is a new point. vertices_determining_current_stack[current_ps.depth] = i new_vertex = 1 @@ -1106,9 +521,8 @@ subgroup_primary_orbit_size += 1 # Look for a new point to split. i += 1 - while i < n and not bitset_check(vertices_to_split[current_ps.depth], i): - i += 1 - if i < n: + i = bitset_next(vertices_to_split[current_ps.depth], i) + if i != -1: # There is a new point. vertices_determining_current_stack[current_ps.depth] = i if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i: @@ -1139,7 +553,7 @@ current_kids_are_same = current_ps.depth + 1 if first_and_current_indicator_same > current_ps.depth: first_and_current_indicator_same = current_ps.depth - if canonical_label and label_and_current_indicator_same > current_ps.depth: + if canonical_label and label_and_current_indicator_same >= current_ps.depth: label_and_current_indicator_same = current_ps.depth compared_current_and_label_indicators = 0 @@ -1149,6 +563,7 @@ while 1: i = current_ps.depth current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by) + PS_move_all_mins_to_front(current_ps) if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]: first_and_current_indicator_same += 1 if canonical_label: @@ -1177,7 +592,7 @@ if discrete: if current_ps.depth == first_and_current_indicator_same: PS_get_perm_from(current_ps, first_ps, permutation) - if compare_structures(permutation, id_perm, S) == 0: + if compare_structures(permutation, id_perm, S, S) == 0: automorphism = 1 if not automorphism: if (not canonical_label) or (compared_current_and_label_indicators < 0): @@ -1188,7 +603,7 @@ if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth): update_label = 1 else: - i = compare_structures(current_ps.entries, label_ps.entries, S) + i = compare_structures(current_ps.entries, label_ps.entries, S, S) if i > 0: update_label = 1 elif i < 0: diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/data_structures_pxd.pxi Mon Sep 15 17:34:52 2008 -0700 @@ -0,0 +1,33 @@ + +#***************************************************************************** +# Copyright (C) 2006 - 2008 Robert L. Miller +# +# Distributed under the terms of the GNU General Public License (GPL) +# http://www.gnu.org/licenses/ +#***************************************************************************** + +include '../../../ext/cdefs.pxi' +include '../../../ext/stdsage.pxi' +include '../../../misc/bitset_pxd.pxi' + +cdef struct OrbitPartition: + # Disjoint set data structure for representing the orbits of the generators + # so far found. Also keeps track of the minimum elements of the cells and + # their sizes. + int degree + int *parent + int *rank + int *mcr # minimum cell representatives - only valid at the root of a cell + int *size # also only valid at the root of a cell + +cdef struct PartitionStack: + # Representation of a node of the search tree. A sequence of partitions of + # length depth + 1, each of which is finer than the last. Partition k is + # represented as PS.entries in order, broken up immediately after each + # entry of levels which is at most k. + int *entries + int *levels + int depth + int degree + + diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/data_structures_pyx.pxi Mon Sep 15 17:34:52 2008 -0700 @@ -0,0 +1,619 @@ +r""" +Data structures + +This module implements TODO... + +REFERENCE: + + [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium, + Vol. 30 (1981), pp. 45-87. + +""" + +#***************************************************************************** +# Copyright (C) 2006 - 2008 Robert L. Miller +# +# Distributed under the terms of the GNU General Public License (GPL) +# http://www.gnu.org/licenses/ +#***************************************************************************** + +include '../../../misc/bitset.pxi' + +cdef inline int min(int a, int b): + if a < b: return a + else: return b + +cdef inline int max(int a, int b): + if a > b: return a + else: return b + +# OrbitPartitions + +cdef inline OrbitPartition *OP_new(int n): + """ + Allocate and return a pointer to a new OrbitPartition of degree n. Returns a + null pointer in the case of an allocation failure. + """ + cdef int i + cdef OrbitPartition *OP = \ + sage_malloc(sizeof(OrbitPartition)) + if OP is NULL: + return OP + OP.degree = n + OP.parent = sage_malloc( n * sizeof(int) ) + OP.rank = sage_malloc( n * sizeof(int) ) + OP.mcr = sage_malloc( n * sizeof(int) ) + OP.size = sage_malloc( n * sizeof(int) ) + if OP.parent is NULL or OP.rank is NULL or \ + OP.mcr is NULL or OP.size is NULL: + if OP.parent is not NULL: sage_free(OP.parent) + if OP.rank is not NULL: sage_free(OP.rank) + if OP.mcr is not NULL: sage_free(OP.mcr) + if OP.size is not NULL: sage_free(OP.size) + return NULL + for i from 0 <= i < n: + OP.parent[i] = i + OP.rank[i] = 0 + OP.mcr[i] = i + OP.size[i] = 1 + return OP + +cdef inline int OP_dealloc(OrbitPartition *OP): + sage_free(OP.parent) + sage_free(OP.rank) + sage_free(OP.mcr) + sage_free(OP.size) + sage_free(OP) + return 0 + +cdef inline int OP_find(OrbitPartition *OP, int n): + """ + Report the representative ("root") of the cell which contains n. + """ + if OP.parent[n] == n: + return n + else: + OP.parent[n] = OP_find(OP, OP.parent[n]) + return OP.parent[n] + +cdef inline int OP_join(OrbitPartition *OP, int m, int n): + """ + Join the cells containing m and n, if they are different. + """ + cdef int m_root = OP_find(OP, m) + cdef int n_root = OP_find(OP, n) + if OP.rank[m_root] > OP.rank[n_root]: + OP.parent[n_root] = m_root + OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root]) + OP.size[m_root] += OP.size[n_root] + elif OP.rank[m_root] < OP.rank[n_root]: + OP.parent[m_root] = n_root + OP.mcr[n_root] = min(OP.mcr[m_root], OP.mcr[n_root]) + OP.size[n_root] += OP.size[m_root] + elif m_root != n_root: + OP.parent[n_root] = m_root + OP.mcr[m_root] = min(OP.mcr[m_root], OP.mcr[n_root]) + OP.size[m_root] += OP.size[n_root] + OP.rank[m_root] += 1 + +cdef inline int OP_merge_list_perm(OrbitPartition *OP, int *gamma): + """ + Joins the cells of OP which intersect the same orbit of gamma. + + INPUT: + gamma - an integer array representing i -> gamma[i]. + + OUTPUT: + 1 - something changed + 0 - orbits of gamma all contained in cells of OP + """ + cdef int i, i_root, gamma_i_root, changed = 0 + for i from 0 <= i < OP.degree: + if gamma[i] == i: continue + i_root = OP_find(OP, i) + gamma_i_root = OP_find(OP, gamma[i]) + if i_root != gamma_i_root: + changed = 1 + OP_join(OP, i_root, gamma_i_root) + return changed + +def OP_represent(int n=9, merges=[(0,1),(2,3),(3,4)], perm=[1,2,0,4,3,6,7,5,8]): + """ + Demonstration and testing. + + DOCTEST: + sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \ + ... import OP_represent + sage: OP_represent() + Allocating OrbitPartition... + Allocation passed. + Checking that each element reports itself as its root. + Each element reports itself as its root. + Merging: + Merged 0 and 1. + Merged 2 and 3. + Merged 3 and 4. + Done merging. + Finding: + 0 -> 0, root: size=2, mcr=0, rank=1 + 1 -> 0 + 2 -> 2, root: size=3, mcr=2, rank=1 + 3 -> 2 + 4 -> 2 + 5 -> 5, root: size=1, mcr=5, rank=0 + 6 -> 6, root: size=1, mcr=6, rank=0 + 7 -> 7, root: size=1, mcr=7, rank=0 + 8 -> 8, root: size=1, mcr=8, rank=0 + Allocating array to test merge_perm. + Allocation passed. + Merging permutation: [1, 2, 0, 4, 3, 6, 7, 5, 8] + Done merging. + Finding: + 0 -> 0, root: size=5, mcr=0, rank=2 + 1 -> 0 + 2 -> 0 + 3 -> 0 + 4 -> 0 + 5 -> 5, root: size=3, mcr=5, rank=1 + 6 -> 5 + 7 -> 5 + 8 -> 8, root: size=1, mcr=8, rank=0 + Deallocating OrbitPartition. + Done. + + """ + cdef int i + print "Allocating OrbitPartition..." + cdef OrbitPartition *OP = OP_new(n) + if OP is NULL: + print "Allocation failed!" + return + print "Allocation passed." + print "Checking that each element reports itself as its root." + good = True + for i from 0 <= i < n: + if not OP_find(OP, i) == i: + print "Failed at i = %d!"%i + good = False + if good: print "Each element reports itself as its root." + print "Merging:" + for i,j in merges: + OP_join(OP, i, j) + print "Merged %d and %d."%(i,j) + print "Done merging." + print "Finding:" + for i from 0 <= i < n: + j = OP_find(OP, i) + s = "%d -> %d"%(i, j) + if i == j: + s += ", root: size=%d, mcr=%d, rank=%d"%\ + (OP.size[i], OP.mcr[i], OP.rank[i]) + print s + print "Allocating array to test merge_perm." + cdef int *gamma = sage_malloc( n * sizeof(int) ) + if gamma is NULL: + print "Allocation failed!" + OP_dealloc(OP) + return + print "Allocation passed." + for i from 0 <= i < n: + gamma[i] = perm[i] + print "Merging permutation: %s"%perm + OP_merge_list_perm(OP, gamma) + print "Done merging." + print "Finding:" + for i from 0 <= i < n: + j = OP_find(OP, i) + s = "%d -> %d"%(i, j) + if i == j: + s += ", root: size=%d, mcr=%d, rank=%d"%\ + (OP.size[i], OP.mcr[i], OP.rank[i]) + print s + print "Deallocating OrbitPartition." + sage_free(gamma) + OP_dealloc(OP) + print "Done." + +# PartitionStacks + +cdef inline PartitionStack *PS_new(int n, bint unit_partition): + """ + Allocate and return a pointer to a new PartitionStack of degree n. Returns a + null pointer in the case of an allocation failure. + """ + cdef int i + cdef PartitionStack *PS = \ + sage_malloc(sizeof(PartitionStack)) + if PS is NULL: + return PS + PS.entries = sage_malloc( n * sizeof(int) ) + PS.levels = sage_malloc( n * sizeof(int) ) + if PS.entries is NULL or PS.levels is NULL: + if PS.entries is not NULL: sage_free(PS.entries) + if PS.levels is not NULL: sage_free(PS.levels) + return NULL + PS.depth = 0 + PS.degree = n + if unit_partition: + for i from 0 <= i < n-1: + PS.entries[i] = i + PS.levels[i] = n + PS.entries[n-1] = n-1 + PS.levels[n-1] = -1 + return PS + +cdef inline PartitionStack *PS_copy(PartitionStack *PS): + """ + Allocate and return a pointer to a copy of PartitionStack PS. Returns a null + pointer in the case of an allocation failure. + """ + cdef int i + cdef PartitionStack *PS2 = \ + sage_malloc(sizeof(PartitionStack)) + if PS2 is NULL: + return PS2 + PS2.entries = sage_malloc( PS.degree * sizeof(int) ) + PS2.levels = sage_malloc( PS.degree * sizeof(int) ) + if PS2.entries is NULL or PS2.levels is NULL: + if PS2.entries is not NULL: sage_free(PS2.entries) + if PS2.levels is not NULL: sage_free(PS2.levels) + return NULL + PS_copy_from_to(PS, PS2) + return PS2 + +cdef inline int PS_copy_from_to(PartitionStack *PS, PartitionStack *PS2): + """ + Copy all data from PS to PS2. + """ + PS2.depth = PS.depth + PS2.degree = PS.degree + for i from 0 <= i < PS.degree: + PS2.entries[i] = PS.entries[i] + PS2.levels[i] = PS.levels[i] + +cdef inline PartitionStack *PS_from_list(object L): + """ + Allocate and return a pointer to a PartitionStack representing L. Returns a + null pointer in the case of an allocation failure. + """ + cdef int cell, i, num_cells = len(L), cur_start = 0, cur_len, n = 0 + for cell from 0 <= cell < num_cells: + n += len(L[cell]) + cdef PartitionStack *PS = PS_new(n, 0) + cell = 0 + if PS is NULL: + return PS + while 1: + cur_len = len(L[cell]) + for i from 0 <= i < cur_len: + PS.entries[cur_start + i] = L[cell][i] + PS.levels[cur_start + i] = n + PS_move_min_to_front(PS, cur_start, cur_start+cur_len-1) + cur_start += cur_len + cell += 1 + if cell == num_cells: + PS.levels[cur_start-1] = -1 + break + PS.levels[cur_start-1] = 0 + return PS + +cdef inline int PS_dealloc(PartitionStack *PS): + sage_free(PS.entries) + sage_free(PS.levels) + sage_free(PS) + return 0 + +cdef inline int PS_print(PartitionStack *PS): + """ + Print a visual representation of PS. + """ + cdef int i + for i from 0 <= i < PS.depth + 1: + PS_print_partition(PS, i) + print + +cdef inline int PS_print_partition(PartitionStack *PS, int k): + """ + Print the partition at depth k. + """ + s = '(' + for i from 0 <= i < PS.degree: + s += str(PS.entries[i]) + if PS.levels[i] <= k: + s += '|' + else: + s += ',' + s = s[:-1] + ')' + print s + +cdef inline bint PS_is_discrete(PartitionStack *PS): + """ + Returns whether the deepest partition consists only of singleton cells. + """ + cdef int i + for i from 0 <= i < PS.degree: + if PS.levels[i] > PS.depth: + return 0 + return 1 + +cdef inline int PS_num_cells(PartitionStack *PS): + """ + Returns the number of cells. + """ + cdef int i, ncells = 0 + for i from 0 <= i < PS.degree: + if PS.levels[i] <= PS.depth: + ncells += 1 + return ncells + +cdef inline int PS_move_min_to_front(PartitionStack *PS, int start, int end): + """ + Makes sure that the first element of the segment of entries i with + start <= i <= end is minimal. + """ + cdef int i, min_loc = start, minimum = PS.entries[start] + for i from start < i <= end: + if PS.entries[i] < minimum: + min_loc = i + minimum = PS.entries[i] + if min_loc != start: + PS.entries[min_loc] = PS.entries[start] + PS.entries[start] = minimum + +cdef inline bint PS_is_mcr(PartitionStack *PS, int m): + """ + Returns whether PS.elements[m] (not m!) is the smallest element of its cell. + """ + return m == 0 or PS.levels[m-1] <= PS.depth + +cdef inline bint PS_is_fixed(PartitionStack *PS, int m): + """ + Returns whether PS.elements[m] (not m!) is in a singleton cell, assuming + PS_is_mcr(PS, m) is already true. + """ + return PS.levels[m] <= PS.depth + +cdef inline int PS_clear(PartitionStack *PS): + """ + Sets the current partition to the first shallower one, i.e. forgets about + boundaries between cells that are new to the current level. + """ + cdef int i, cur_start = 0 + for i from 0 <= i < PS.degree: + if PS.levels[i] >= PS.depth: + PS.levels[i] += 1 + if PS.levels[i] < PS.depth: + PS_move_min_to_front(PS, cur_start, i) + cur_start = i+1 + +cdef inline int PS_move_all_mins_to_front(PartitionStack *PS): + """ + Move minimal cell elements to the front of each cell. + """ + cdef int i, cur_start = 0 + for i from 0 <= i < PS.degree: + if PS.levels[i] <= PS.depth: + PS_move_min_to_front(PS, cur_start, i) + cur_start = i+1 + +cdef inline int PS_split_point(PartitionStack *PS, int v): + """ + Detaches the point v from the cell it is in, putting the singleton cell + of just v in front. Returns the position where v is now located. + """ + cdef int i = 0, index_of_v + while PS.entries[i] != v: + i += 1 + index_of_v = i + while PS.levels[i] > PS.depth: + i += 1 + if (index_of_v == 0 or PS.levels[index_of_v-1] <= PS.depth) \ + and PS.levels[index_of_v] > PS.depth: + # if v is first (make sure v is not already alone) + PS_move_min_to_front(PS, index_of_v+1, i) + PS.levels[index_of_v] = PS.depth + return index_of_v + else: + # If v is not at front, i.e. v is not minimal in its cell, + # then move_min_to_front is not necessary since v will swap + # with the first before making its own cell, leaving it at + # the front of the other. + i = index_of_v + while i != 0 and PS.levels[i-1] > PS.depth: + i -= 1 + PS.entries[index_of_v] = PS.entries[i+1] # move the second element to v + PS.entries[i+1] = PS.entries[i] # move the first (min) to second + PS.entries[i] = v # place v first + PS.levels[i] = PS.depth + return i + +cdef inline int PS_first_smallest(PartitionStack *PS, bitset_t b): + """ + Find the first occurrence of the smallest cell of size greater than one, + store its entries to b, and return its minimum element. + """ + cdef int i = 0, j = 0, location = 0, n = PS.degree + bitset_zero(b) + while 1: + if PS.levels[i] <= PS.depth: + if i != j and n > i - j + 1: + n = i - j + 1 + location = j + j = i + 1 + if PS.levels[i] == -1: break + i += 1 + # location now points to the beginning of the first, smallest, + # nontrivial cell + i = location + while 1: + bitset_flip(b, PS.entries[i]) + if PS.levels[i] <= PS.depth: break + i += 1 + return PS.entries[location] + +cdef inline int PS_get_perm_from(PartitionStack *PS1, PartitionStack *PS2, int *gamma): + """ + Store the permutation determined by PS2[i] -> PS1[i] for each i, where PS[i] + denotes the entry of the ith cell of the discrete partition PS. + """ + cdef int i + for i from 0 <= i < PS1.degree: + gamma[PS2.entries[i]] = PS1.entries[i] + +def PS_represent(partition=[[6],[3,4,8,7],[1,9,5],[0,2]], splits=[6,1,8,2]): + """ + Demonstration and testing. + + DOCTEST: + sage: from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label \ + ... import PS_represent + sage: PS_represent() + Allocating PartitionStack... + Allocation passed: + (0,1,2,3,4,5,6,7,8,9) + Checking that entries are in order and correct level. + Everything seems in order, deallocating. + Deallocated. + Creating PartitionStack from partition [[6], [3, 4, 8, 7], [1, 9, 5], [0, 2]]. + PartitionStack's data: + entries -> [6, 3, 4, 8, 7, 1, 9, 5, 0, 2] + levels -> [0, 10, 10, 10, 0, 10, 10, 0, 10, -1] + depth = 0, degree = 10 + (6|3,4,8,7|1,9,5|0,2) + Checking PS_is_discrete: + False + Checking PS_num_cells: + 4 + Checking PS_is_mcr, min cell reps are: + [6, 3, 1, 0] + Checking PS_is_fixed, fixed elements are: + [6] + Copying PartitionStack: + (6|3,4,8,7|1,9,5|0,2) + Checking for consistency. + Everything is consistent. + Clearing copy: + (0,3,4,8,7,1,9,5,6,2) + Splitting point 6 from original: + 0 + (6|3,4,8,7|1,9,5|0,2) + Splitting point 1 from original: + 5 + (6|3,4,8,7|1|5,9|0,2) + Splitting point 8 from original: + 1 + (6|8|3,4,7|1|5,9|0,2) + Splitting point 2 from original: + 8 + (6|8|3,4,7|1|5,9|2|0) + Getting permutation from PS2->PS: + [6, 1, 0, 8, 3, 9, 2, 7, 4, 5] + Finding first smallest: + Minimal element is 5, bitset is: + 0000010001 + Deallocating PartitionStacks. + Done. + + """ + cdef int i, n = sum([len(cell) for cell in partition]), *gamma + cdef bitset_t b + print "Allocating PartitionStack..." + cdef PartitionStack *PS = PS_new(n, 1), *PS2 + if PS is NULL: + print "Allocation failed!" + return + print "Allocation passed:" + PS_print(PS) + print "Checking that entries are in order and correct level." + good = True + for i from 0 <= i < n-1: + if not (PS.entries[i] == i and PS.levels[i] == n): + print "Failed at i = %d!"%i + print PS.entries[i], PS.levels[i], i, n + good = False + if not (PS.entries[n-1] == n-1 and PS.levels[n-1] == -1): + print "Failed at i = %d!"%(n-1) + good = False + if not PS.degree == n or not PS.depth == 0: + print "Incorrect degree or depth!" + good = False + if good: print "Everything seems in order, deallocating." + PS_dealloc(PS) + print "Deallocated." + print "Creating PartitionStack from partition %s."%partition + PS = PS_from_list(partition) + print "PartitionStack's data:" + print "entries -> %s"%[PS.entries[i] for i from 0 <= i < n] + print "levels -> %s"%[PS.levels[i] for i from 0 <= i < n] + print "depth = %d, degree = %d"%(PS.depth,PS.degree) + PS_print(PS) + print "Checking PS_is_discrete:" + print "True" if PS_is_discrete(PS) else "False" + print "Checking PS_num_cells:" + print PS_num_cells(PS) + print "Checking PS_is_mcr, min cell reps are:" + L = [PS.entries[i] for i from 0 <= i < n if PS_is_mcr(PS, i)] + print L + print "Checking PS_is_fixed, fixed elements are:" + print [PS.entries[l] for l in L if PS_is_fixed(PS, l)] + print "Copying PartitionStack:" + PS2 = PS_copy(PS) + PS_print(PS2) + print "Checking for consistency." + good = True + for i from 0 <= i < n: + if PS.entries[i] != PS2.entries[i] or PS.levels[i] != PS2.levels[i]: + print "Failed at i = %d!"%i + good = False + if PS.degree != PS2.degree or PS.depth != PS2.depth: + print "Failure with degree or depth!" + good = False + if good: + print "Everything is consistent." + print "Clearing copy:" + PS_clear(PS2) + PS_print(PS2) + for s in splits: + print "Splitting point %d from original:"%s + print PS_split_point(PS, s) + PS_print(PS) + print "Getting permutation from PS2->PS:" + gamma = sage_malloc(n * sizeof(int)) + PS_get_perm_from(PS, PS2, gamma) + print [gamma[i] for i from 0 <= i < n] + sage_free(gamma) + print "Finding first smallest:" + bitset_init(b, n) + i = PS_first_smallest(PS, b) + print "Minimal element is %d, bitset is:"%i + print bitset_string(b) + bitset_clear(b) + print "Deallocating PartitionStacks." + PS_dealloc(PS) + PS_dealloc(PS2) + print "Done." + +# Functions + +cdef inline int split_point_and_refine(PartitionStack *PS, int v, object S, + int (*refine_and_return_invariant)\ + (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len), + int *cells_to_refine_by): + """ + Make the partition stack one longer by copying the last partition in the + stack, split off a given point, and refine. Return the invariant given by + the refinement function. + + INPUT: + PS -- the partition stack to refine + v -- the point to split + S -- the structure + refine_and_return_invariant -- the refinement function provided + cells_to_refine_by -- an array, contents ignored + + """ + PS.depth += 1 + PS_clear(PS) + cells_to_refine_by[0] = PS_split_point(PS, v) + return refine_and_return_invariant(PS, S, cells_to_refine_by, 1) + diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/double_coset.pxd --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/double_coset.pxd Mon Sep 15 17:34:52 2008 -0700 @@ -0,0 +1,21 @@ + +#***************************************************************************** +# Copyright (C) 2006 - 2008 Robert L. Miller +# +# Distributed under the terms of the GNU General Public License (GPL) +# http://www.gnu.org/licenses/ +#***************************************************************************** + +include '../../../ext/cdefs.pxi' +include '../../../ext/stdsage.pxi' +include 'data_structures_pxd.pxi' # includes bitsets + +from sage.rings.integer cimport Integer + +cdef int *double_coset( object, object, int **, int *, int, + bint (*)(PartitionStack *, object), + int (*)(PartitionStack *, object, int *, int), + int (*)(int *, int *, object, object)) + + + diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/double_coset.pyx --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/double_coset.pyx Mon Sep 15 17:34:52 2008 -0700 @@ -0,0 +1,595 @@ +r""" +Double cosets + +This module implements a general algorithm for computing double coset problems +for pairs of objects. The class of objects in question must be some kind +of structure for which an isomorphism is a permutation in $S_n$ for some $n$, +which we call here the order of the object. Given objects $X$ and $Y$, +the program returns an isomorphism in list permutation form if $X \cong Y$, and +a NULL pointer otherwise. + +In order to take advantage of the algorithms in this module for a specific kind +of object, one must implement (in Cython) three functions which will be specific +to the kind of objects in question. Pointers to these functions are passed to +the main function of the module, which is \code{double_coset}. For specific +examples of implementations of these functions, see any of the files in +\code{sage.groups.perm_gps.partn_ref} beginning with "refinement." They are: + +A. \code{refine_and_return_invariant}: + + Signature: + + \code{int refine_and_return_invariant(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len)} + + This function should split up cells in the partition at the top of the + partition stack in such a way that any automorphism that respects the + partition also respects the resulting partition. The array + cells_to_refine_by is a list of the beginning positions of some cells which + have been changed since the last refinement. It is not necessary to use + this in an implementation of this function, but it will affect performance. + One should consult \code{refinement_graphs} for more details and ideas for + particular implementations. + + Output: + + An integer $I$ invariant under the orbits of $S_n$. That is, if + $\gamma \in S_n$, then + $$ I(G, PS, cells_to_refine_by) = I( \gamma(G), \gamma(PS), \gamma(cells_to_refine_by) ) .$$ + + +B. \code{compare_structures}: + + Signature: + + \code{int compare_structures(int *gamma_1, int *gamma_2, object S)} + + This function must implement a total ordering on the set of objects of fixed + order. Return -1 if \code{gamma_1(S) < gamma_2(S)}, 0 if + \code{gamma_1(S) == gamma_2(S)}, 1 if \code{gamma_1(S) > gamma_2(S)}. + +C. \code{all_children_are_equivalent}: + + Signature: + + \code{bint all_children_are_equivalent(PartitionStack *PS, object S)} + + This function must return False unless it is the case that each discrete + partition finer than the top of the partition stack is equivalent to the + others under some automorphism of S. The converse need not hold: if this is + indeed the case, it still may return False. This function is originally used + as a consequence of Lemma 2.25 in [1]. + +DOCTEST: + sage: import sage.groups.perm_gps.partn_ref.double_coset + +REFERENCE: + + [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium, + Vol. 30 (1981), pp. 45-87. + + [2] Leon, Jeffrey. Permutation Group Algorithms Based on Partitions, I: + Theory and Algorithms. J. Symbolic Computation, Vol. 12 (1991), pp. + 533-583. + +""" + +#***************************************************************************** +# Copyright (C) 2006 - 2008 Robert L. Miller +# +# Distributed under the terms of the GNU General Public License (GPL) +# http://www.gnu.org/licenses/ +#***************************************************************************** + +include 'data_structures_pyx.pxi' # includes bitsets + +cdef inline int cmp(int a, int b): + if a < b: return -1 + elif a == b: return 0 + else: return 1 + +cdef inline bint stacks_are_equivalent(PartitionStack *PS1, PartitionStack *PS2): + cdef int i, j, depth = min(PS1.depth, PS2.depth) + for i from 0 <= i < PS1.degree: + j = cmp(PS1.levels[i], PS2.levels[i]) + if j == 0: continue + if ( (j < 0 and PS1.levels[i] <= depth and PS2.levels[i] > depth) + or (j > 0 and PS2.levels[i] <= depth and PS1.levels[i] > depth) ): + return 0 + return 1 + +# Functions + +cdef int *double_coset(object S1, object S2, int **partition1, int *ordering2, + int n, bint (*all_children_are_equivalent)(PartitionStack *PS, object S), + int (*refine_and_return_invariant)\ + (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len), + int (*compare_structures)(int *gamma_1, int *gamma_2, object S1, object S2) ): + """ + Traverse the search space for double coset calculation. + + INPUT: + S1, S2 -- pointers to the structures + partition1 -- array representing a partition of the points of S1 + ordering2 -- an ordering of the points of S2 representing a second partition + n -- the number of points (points are assumed to be 0,1,...,n-1) + all_children_are_equivalent -- pointer to a function + INPUT: + PS -- pointer to a partition stack + S -- pointer to the structure + OUTPUT: + bint -- returns True if it can be determined that all refinements below + the current one will result in an equivalent discrete partition + refine_and_return_invariant -- pointer to a function + INPUT: + PS -- pointer to a partition stack + S -- pointer to the structure + alpha -- an array consisting of numbers, which indicate the starting + positions of the cells to refine against (will likely be modified) + OUTPUT: + int -- returns an invariant under application of arbitrary permutations + compare_structures -- pointer to a function + INPUT: + gamma_1, gamma_2 -- (list) permutations of the points of S1 and S2 + S1, S2 -- pointers to the structures + OUTPUT: + int -- 0 if gamma_1(S1) = gamma_2(S2), otherwise -1 or 1 (see docs for cmp), + such that the set of all structures is well-ordered + OUTPUT: + If S1 and S2 are isomorphic, a pointer to an integer array representing an + isomorphism. Otherwise, a NULL pointer. + + """ + cdef PartitionStack *current_ps, *first_ps, *left_ps + cdef int first_meets_current = -1 + cdef int current_kids_are_same = 1 + cdef int first_kids_are_same + + cdef int *indicators + + cdef OrbitPartition *orbits_of_subgroup + cdef int subgroup_primary_orbit_size = 0 + cdef int minimal_in_primary_orbit + + cdef bitset_t *fixed_points_of_generators # i.e. fp + cdef bitset_t *minimal_cell_reps_of_generators # i.e. mcr + cdef int len_of_fp_and_mcr = 100 + cdef int index_in_fp_and_mcr = -1 + + cdef bitset_t *vertices_to_split + cdef bitset_t vertices_have_been_reduced + cdef int *permutation, *id_perm, *cells_to_refine_by + cdef int *vertices_determining_current_stack + + cdef int *output + + cdef int i, j, k + cdef bint discrete, automorphism, update_label + cdef bint backtrack, new_vertex, narrow, mem_err = 0 + + if n == 0: + return NULL + + indicators = sage_malloc(n * sizeof(int)) + + fixed_points_of_generators = sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) ) + minimal_cell_reps_of_generators = sage_malloc( len_of_fp_and_mcr * sizeof(bitset_t) ) + + vertices_to_split = sage_malloc( n * sizeof(bitset_t) ) + permutation = sage_malloc( n * sizeof(int) ) + id_perm = sage_malloc( n * sizeof(int) ) + cells_to_refine_by = sage_malloc( n * sizeof(int) ) + vertices_determining_current_stack = sage_malloc( n * sizeof(int) ) + + output = sage_malloc( n * sizeof(int) ) + + current_ps = PS_new(n, 0) + left_ps = PS_new(n, 0) + orbits_of_subgroup = OP_new(n) + + # Check for allocation failures: + if indicators is NULL or \ + fixed_points_of_generators is NULL or \ + minimal_cell_reps_of_generators is NULL or \ + vertices_to_split is NULL or \ + permutation is NULL or \ + id_perm is NULL or \ + cells_to_refine_by is NULL or \ + vertices_determining_current_stack is NULL or \ + current_ps is NULL or \ + orbits_of_subgroup is NULL or \ + output is NULL: + if indicators is not NULL: + sage_free(indicators) + if fixed_points_of_generators is not NULL: + sage_free(fixed_points_of_generators) + if minimal_cell_reps_of_generators is not NULL: + sage_free(minimal_cell_reps_of_generators) + if vertices_to_split is not NULL: + sage_free(vertices_to_split) + if permutation is not NULL: + sage_free(permutation) + if id_perm is not NULL: + sage_free(id_perm) + if cells_to_refine_by is not NULL: + sage_free(cells_to_refine_by) + if vertices_determining_current_stack is not NULL: + sage_free(vertices_determining_current_stack) + if current_ps is not NULL: + PS_dealloc(current_ps) + if orbits_of_subgroup is not NULL: + OP_dealloc(orbits_of_subgroup) + if output is not NULL: + sage_free(output) + raise MemoryError + + # Initialize bitsets, checking for allocation failures: + cdef bint succeeded = 1 + for i from 0 <= i < len_of_fp_and_mcr: + try: + bitset_init(fixed_points_of_generators[i], n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < i: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + break + try: + bitset_init(minimal_cell_reps_of_generators[i], n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < i: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + bitset_clear(fixed_points_of_generators[i]) + break + if succeeded: + for i from 0 <= i < n: + try: + bitset_init(vertices_to_split[i], n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < i: + bitset_clear(vertices_to_split[j]) + for j from 0 <= j < len_of_fp_and_mcr: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + break + if succeeded: + try: + bitset_init(vertices_have_been_reduced, n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < n: + bitset_clear(vertices_to_split[j]) + for j from 0 <= j < len_of_fp_and_mcr: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + if not succeeded: + sage_free(indicators) + sage_free(fixed_points_of_generators) + sage_free(minimal_cell_reps_of_generators) + sage_free(vertices_to_split) + sage_free(permutation) + sage_free(id_perm) + sage_free(cells_to_refine_by) + sage_free(vertices_determining_current_stack) + PS_dealloc(current_ps) + PS_dealloc(left_ps) + OP_dealloc(orbits_of_subgroup) + raise MemoryError + + bitset_zero(vertices_have_been_reduced) + + # set up the identity permutation + for i from 0 <= i < n: + id_perm[i] = i + if ordering2 is NULL: + ordering2 = id_perm + + # Copy data of partition to left_ps, and + # ordering of that partition to current_ps. + i = 0 + j = 0 + while partition1[i] is not NULL: + k = 0 + while partition1[i][k] != -1: + left_ps.entries[j+k] = partition1[i][k] + left_ps.levels[j+k] = n + current_ps.entries[j+k] = ordering2[j+k] + current_ps.levels[j+k] = n + k += 1 + left_ps.levels[j+k-1] = 0 + current_ps.levels[j+k-1] = 0 + PS_move_min_to_front(current_ps, j, j+k-1) + j += k + i += 1 + left_ps.levels[j-1] = -1 + current_ps.levels[j-1] = -1 + left_ps.depth = 0 + current_ps.depth = 0 + left_ps.degree = n + current_ps.degree = n + + # default values of "infinity" + for i from 0 <= i < n: + indicators[i] = -1 + + cdef bint possible = 1 + cdef bint unknown = 1 + + # Our first refinement needs to check every cell of the partition, + # so cells_to_refine_by needs to be a list of pointers to each cell. + j = 1 + cells_to_refine_by[0] = 0 + for i from 0 < i < n: + if left_ps.levels[i-1] == 0: + cells_to_refine_by[j] = i + j += 1 + + k = refine_and_return_invariant(left_ps, S1, cells_to_refine_by, j) + + j = 1 + cells_to_refine_by[0] = 0 + for i from 0 < i < n: + if current_ps.levels[i-1] == 0: + cells_to_refine_by[j] = i + j += 1 + j = refine_and_return_invariant(current_ps, S2, cells_to_refine_by, j) + if k != j: + possible = 0; unknown = 0 + elif not stacks_are_equivalent(left_ps, current_ps): + possible = 0; unknown = 0 + else: + PS_move_all_mins_to_front(current_ps) + + first_ps = NULL + # Refine down to a discrete partition + while not PS_is_discrete(left_ps) and possible: + k = PS_first_smallest(left_ps, vertices_to_split[left_ps.depth]) + i = left_ps.depth + indicators[i] = split_point_and_refine(left_ps, k, S1, refine_and_return_invariant, cells_to_refine_by) + vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth]) + bitset_unset(vertices_have_been_reduced, current_ps.depth) + while 1: + j = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by) + if indicators[i] != j: + possible = 0 + elif not stacks_are_equivalent(left_ps, current_ps): + possible = 0 + else: + PS_move_all_mins_to_front(current_ps) + if not possible: + j = vertices_determining_current_stack[i] + 1 + j = bitset_next(vertices_to_split[i], j) + if j == -1: + break + else: + possible = 1 + vertices_determining_current_stack[i] = j + current_ps.depth -= 1 # reset for next refinement + else: break + if not all_children_are_equivalent(current_ps, S2): + current_kids_are_same = current_ps.depth + 1 + + if possible: + if compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0: + unknown = 0 + else: + first_meets_current = current_ps.depth + first_kids_are_same = current_ps.depth + first_ps = PS_copy(current_ps) + current_ps.depth -= 1 + + # Main loop: + while possible and unknown and current_ps.depth != -1: + + # I. Search for a new vertex to split, and update subgroup information + new_vertex = 0 + if current_ps.depth > first_meets_current: + # If we are not at a node of the first stack, reduce size of + # vertices_to_split by using the symmetries we already know. + if not bitset_check(vertices_have_been_reduced, current_ps.depth): + for i from 0 <= i <= index_in_fp_and_mcr: + j = 0 + while j < current_ps.depth and bitset_check(fixed_points_of_generators[i], vertices_determining_current_stack[j]): + j += 1 + # If each vertex split so far is fixed by generator i, + # then remove elements of vertices_to_split which are + # not minimal in their orbits under generator i. + if j == current_ps.depth: + for k from 0 <= k < n: + if bitset_check(vertices_to_split[current_ps.depth], k) and not bitset_check(minimal_cell_reps_of_generators[i], k): + bitset_flip(vertices_to_split[current_ps.depth], k) + bitset_flip(vertices_have_been_reduced, current_ps.depth) + # Look for a new point to split. + i = vertices_determining_current_stack[current_ps.depth] + 1 + i = bitset_next(vertices_to_split[current_ps.depth], i) + if i != -1: + # There is a new point. + vertices_determining_current_stack[current_ps.depth] = i + new_vertex = 1 + else: + # No new point: backtrack. + current_ps.depth -= 1 + else: + # If we are at a node of the first stack, the above reduction + # will not help. Also, we must update information about + # primary orbits here. + if current_ps.depth < first_meets_current: + # If we are done searching under this part of the first + # stack, then first_meets_current is one higher, and we + # are looking at a new primary orbit (corresponding to a + # larger subgroup in the stabilizer chain). + first_meets_current = current_ps.depth + for i from 0 <= i < n: + if bitset_check(vertices_to_split[current_ps.depth], i): + minimal_in_primary_orbit = i + break + while 1: + i = vertices_determining_current_stack[current_ps.depth] + # This was the last point to be split here. + # If it is in the same orbit as minimal_in_primary_orbit, + # then count it as an element of the primary orbit. + if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit): + subgroup_primary_orbit_size += 1 + # Look for a new point to split. + i += 1 + i = bitset_next(vertices_to_split[current_ps.depth], i) + if i != -1: + # There is a new point. + vertices_determining_current_stack[current_ps.depth] = i + if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i: + new_vertex = 1 + break + else: + # No new point: backtrack. + # Note that now, we are backtracking up the first stack. + vertices_determining_current_stack[current_ps.depth] = -1 + # If every choice of point to split gave something in the + # (same!) primary orbit, then all children of the first + # stack at this point are equivalent. + j = 0 + for i from 0 <= i < n: + if bitset_check(vertices_to_split[current_ps.depth], i): + j += 1 + if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1: + first_kids_are_same = current_ps.depth + # Backtrack. + subgroup_primary_orbit_size = 0 + current_ps.depth -= 1 + break + if not new_vertex: + continue + + if current_kids_are_same > current_ps.depth + 1: + current_kids_are_same = current_ps.depth + 1 + + # II. Refine down to a discrete partition, or until + # we leave the part of the tree we are interested in + discrete = 0 + while 1: + i = current_ps.depth + while 1: + k = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S2, refine_and_return_invariant, cells_to_refine_by) + PS_move_all_mins_to_front(current_ps) + if indicators[i] != k: + possible = 0 + elif not stacks_are_equivalent(left_ps, current_ps): + possible = 0 + if not possible: + j = vertices_determining_current_stack[i] + 1 + j = bitset_next(vertices_to_split[i], j) + if j == -1: + break + else: + possible = 1 + vertices_determining_current_stack[i] = j + current_ps.depth -= 1 # reset for next refinement + else: break + if not possible: + break + if PS_is_discrete(current_ps): + break + vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth]) + bitset_unset(vertices_have_been_reduced, current_ps.depth) + if not all_children_are_equivalent(current_ps, S2): + current_kids_are_same = current_ps.depth + 1 + + # III. Check for automorphisms and isomorphisms + automorphism = 0 + if possible: + PS_get_perm_from(current_ps, first_ps, permutation) + if compare_structures(permutation, id_perm, S2, S2) == 0: + automorphism = 1 + if not automorphism and possible: + # if we get here, discrete must be true + if current_ps.depth != left_ps.depth: + possible = 0 + elif compare_structures(left_ps.entries, current_ps.entries, S1, S2) == 0: + unknown = 0 + break # main loop + else: + possible = 0 + if automorphism: + if index_in_fp_and_mcr < len_of_fp_and_mcr - 1: + index_in_fp_and_mcr += 1 + bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr]) + bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + for i from 0 <= i < n: + if permutation[i] == i: + bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i) + bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i) + else: + bitset_unset(fixed_points_of_generators[index_in_fp_and_mcr], i) + k = i + j = permutation[i] + while j != i: + if j < k: k = j + j = permutation[j] + if k == i: + bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i) + else: + bitset_unset(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i) + current_ps.depth = first_meets_current + if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser + if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit: + continue # main loop + if bitset_check(vertices_have_been_reduced, current_ps.depth): + bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + continue # main loop + if not possible: + possible = 1 + i = current_ps.depth + current_ps.depth = min(first_kids_are_same-1, current_kids_are_same-1) + if i == current_kids_are_same: + continue # main loop + if index_in_fp_and_mcr < len_of_fp_and_mcr - 1: + index_in_fp_and_mcr += 1 + bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr]) + bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + j = current_ps.depth + current_ps.depth = i # just for mcr and fixed functions... + for i from 0 <= i < n: + if PS_is_mcr(current_ps, i): + bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i) + if PS_is_fixed(current_ps, i): + bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i) + current_ps.depth = j + if bitset_check(vertices_have_been_reduced, current_ps.depth): + bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + + # End of main loop. + + if possible and not unknown: + for i from 0 <= i < n: + output[left_ps.entries[i]] = current_ps.entries[i] + else: + sage_free(output) + output = NULL + + # Deallocate: + for i from 0 <= i < len_of_fp_and_mcr: + bitset_clear(fixed_points_of_generators[i]) + bitset_clear(minimal_cell_reps_of_generators[i]) + for i from 0 <= i < n: + bitset_clear(vertices_to_split[i]) + bitset_clear(vertices_have_been_reduced) + sage_free(indicators) + sage_free(fixed_points_of_generators) + sage_free(minimal_cell_reps_of_generators) + sage_free(vertices_to_split) + sage_free(permutation) + sage_free(id_perm) + sage_free(cells_to_refine_by) + sage_free(vertices_determining_current_stack) + PS_dealloc(current_ps) + PS_dealloc(left_ps) + OP_dealloc(orbits_of_subgroup) + if first_ps is not NULL: PS_dealloc(first_ps) + + return output + + + + + diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_binary.pxd --- a/sage/groups/perm_gps/partn_ref/refinement_binary.pxd Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/refinement_binary.pxd Mon Sep 15 17:34:52 2008 -0700 @@ -8,10 +8,11 @@ include '../../../ext/cdefs.pxi' include '../../../ext/stdsage.pxi' -include '../../../misc/bitset_pxd.pxi' +include 'data_structures_pxd.pxi' # includes bitsets from sage.rings.integer cimport Integer -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 +from sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return +from double_coset cimport double_coset cdef class BinaryCodeStruct: cdef bitset_s *alpha_is_wd # single bitset of length nwords + degree @@ -36,8 +37,8 @@ cdef int ith_word_nonlinear(BinaryCodeStruct, int, bitset_s *) cdef int refine_by_bip_degree(PartitionStack *, object, int *, int) -cdef int compare_linear_codes(int *, int *, object) -cdef int compare_nonlinear_codes(int *, int *, object) +cdef int compare_linear_codes(int *, int *, object, object) +cdef int compare_nonlinear_codes(int *, int *, object, object) cdef bint all_children_are_equivalent(PartitionStack *, object) cdef inline int word_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *) cdef inline int col_degree(PartitionStack *, BinaryCodeStruct, int, int, PartitionStack *) diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_binary.pyx --- a/sage/groups/perm_gps/partn_ref/refinement_binary.pyx Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/refinement_binary.pyx Mon Sep 15 17:34:52 2008 -0700 @@ -22,7 +22,8 @@ # http://www.gnu.org/licenses/ #***************************************************************************** -include '../../../misc/bitset.pxi' +include 'data_structures_pyx.pxi' # includes bitsets + from sage.matrix.matrix import is_Matrix cdef class LinearBinaryCodeStruct(BinaryCodeStruct): @@ -95,7 +96,6 @@ self.output = NULL self.ith_word = &ith_word_linear - self.first_time = 1 def run(self, partition=None): """ @@ -238,6 +238,7 @@ part[i][j] = partition[i][j] part[i][len(partition[i])] = -1 part[len(partition)] = NULL + self.first_time = 1 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) @@ -295,6 +296,58 @@ if self.output is NULL: self.run() return [self.output.relabeling[i] for i from 0 <= i < self.degree] + + def is_isomorphic(self, LinearBinaryCodeStruct other): + """ + Calculate whether self is isomorphic to other. + + EXAMPLES: + sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct + + sage: B = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])) + sage: C = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,1],[1,1,0,1,1,0]])) + sage: B.is_isomorphic(C) + [0, 1, 2, 5, 3, 4] + + """ + cdef int **part, i, j + cdef int *output, *ordering + partition = [range(self.degree)] + part = sage_malloc((len(partition)+1) * sizeof(int *)) + ordering = sage_malloc(self.degree * sizeof(int)) + if part is NULL or ordering is NULL: + if part is not NULL: sage_free(part) + if ordering is not NULL: sage_free(ordering) + raise MemoryError + for i from 0 <= i < len(partition): + part[i] = sage_malloc((len(partition[i])+1) * sizeof(int)) + if part[i] is NULL: + for j from 0 <= j < i: + sage_free(part[j]) + sage_free(part) + raise MemoryError + for j from 0 <= j < len(partition[i]): + part[i][j] = partition[i][j] + part[i][len(partition[i])] = -1 + part[len(partition)] = NULL + for i from 0 <= i < self.degree: + ordering[i] = i + self.first_time = 1 + other.first_time = 1 + + output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes) + + for i from 0 <= i < len(partition): + sage_free(part[i]) + sage_free(part) + sage_free(ordering) + + if output is NULL: + return False + else: + output_py = [output[i] for i from 0 <= i < self.degree] + sage_free(output) + return output_py def __dealloc__(self): cdef int j @@ -402,7 +455,6 @@ self.output = NULL self.ith_word = &ith_word_nonlinear - self.first_time = 1 def __dealloc__(self): cdef int j @@ -477,6 +529,7 @@ part[i][j] = partition[i][j] part[i][len(partition[i])] = -1 part[len(partition)] = NULL + self.first_time = 1 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) @@ -535,6 +588,58 @@ self.run() return [self.output.relabeling[i] for i from 0 <= i < self.degree] + def is_isomorphic(self, NonlinearBinaryCodeStruct other): + """ + Calculate whether self is isomorphic to other. + + EXAMPLES: + sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct + + sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]])) + sage: C = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,0,0,1,1],[1,1,1,1,0,0]])) + sage: B.is_isomorphic(C) + [2, 3, 0, 1, 4, 5] + + """ + cdef int **part, i, j + cdef int *output, *ordering + partition = [range(self.degree)] + part = sage_malloc((len(partition)+1) * sizeof(int *)) + ordering = sage_malloc(self.degree * sizeof(int)) + if part is NULL or ordering is NULL: + if part is not NULL: sage_free(part) + if ordering is not NULL: sage_free(ordering) + raise MemoryError + for i from 0 <= i < len(partition): + part[i] = sage_malloc((len(partition[i])+1) * sizeof(int)) + if part[i] is NULL: + for j from 0 <= j < i: + sage_free(part[j]) + sage_free(part) + raise MemoryError + for j from 0 <= j < len(partition[i]): + part[i][j] = partition[i][j] + part[i][len(partition[i])] = -1 + part[len(partition)] = NULL + for i from 0 <= i < self.degree: + ordering[i] = i + self.first_time = 1 + other.first_time = 1 + + output = double_coset(self, other, part, ordering, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes) + + for i from 0 <= i < len(partition): + sage_free(part[i]) + sage_free(part) + sage_free(ordering) + + if output is NULL: + return False + else: + output_py = [output[i] for i from 0 <= i < self.degree] + sage_free(output) + return output_py + cdef int ith_word_nonlinear(BinaryCodeStruct self, int i, bitset_s *word): cdef NonlinearBinaryCodeStruct NBCS = self bitset_copy(word, &NBCS.words[i]) @@ -606,6 +711,7 @@ if necessary_to_split_cell: invariant += 8 first_largest_subcell = sort_by_function(col_ps, current_cell, col_degrees, col_counts, col_output, BCS.nwords+1) + invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps) invariant += first_largest_subcell against_index = current_cell_against while against_index < ctrb_len: @@ -622,7 +728,6 @@ r += 1 if r >= i: break - invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps) invariant += (i - current_cell) current_cell = i else: @@ -641,6 +746,7 @@ if necessary_to_split_cell: invariant += 64 first_largest_subcell = sort_by_function(word_ps, current_cell, word_degrees, word_counts, word_output, BCS.degree+1) + invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps) invariant += first_largest_subcell against_index = current_cell_against while against_index < ctrb_len: @@ -658,18 +764,17 @@ r += 1 if r >= i: break - invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps) invariant += (i - current_cell) current_cell = i current_cell_against += 1 return invariant -cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S): +cdef int compare_linear_codes(int *gamma_1, int *gamma_2, object S1, object S2): """ - Compare gamma_1(B) and gamma_2(B). + Compare gamma_1(S1) and gamma_2(S2). - Return return -1 if gamma_1(B) < gamma_2(B), 0 if gamma_1(B) == gamma_2(B), - 1 if gamma_1(B) > gamma_2(B). (Just like the python \code{cmp}) function. + Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2), + 1 if gamma_1(S1) > gamma_2(S2). (Just like the python \code{cmp}) function. Abstractly, what this function does is relabel the basis of B by gamma_1 and gamma_2, run a row reduction on each, and verify that the matrices are the @@ -680,29 +785,30 @@ INPUT: gamma_1, gamma_2 -- list permutations (inverse) - S -- a binary code struct object + S1, S2 -- binary code struct objects """ cdef int i, piv_loc_1, piv_loc_2, cur_col, cur_row=0 cdef bint is_pivot_1, is_pivot_2 - cdef LinearBinaryCodeStruct BCS = S - cdef bitset_s *basis_1 = BCS.scratch_bitsets # len = dim - cdef bitset_s *basis_2 = &BCS.scratch_bitsets[BCS.dimension] # len = dim - cdef bitset_s *pivots = &BCS.scratch_bitsets[2*BCS.dimension] # len 1 - cdef bitset_s *temp = &BCS.scratch_bitsets[2*BCS.dimension+1] # len 1 - for i from 0 <= i < BCS.dimension: - bitset_copy(&basis_1[i], &BCS.basis[i]) - bitset_copy(&basis_2[i], &BCS.basis[i]) + cdef LinearBinaryCodeStruct BCS1 = S1 + cdef LinearBinaryCodeStruct BCS2 = S2 + cdef bitset_s *basis_1 = BCS1.scratch_bitsets # len = dim + cdef bitset_s *basis_2 = &BCS1.scratch_bitsets[BCS1.dimension] # len = dim + cdef bitset_s *pivots = &BCS1.scratch_bitsets[2*BCS1.dimension] # len 1 + cdef bitset_s *temp = &BCS1.scratch_bitsets[2*BCS1.dimension+1] # len 1 + for i from 0 <= i < BCS1.dimension: + bitset_copy(&basis_1[i], &BCS1.basis[i]) + bitset_copy(&basis_2[i], &BCS2.basis[i]) bitset_zero(pivots) - for cur_col from 0 <= cur_col < BCS.degree: + for cur_col from 0 <= cur_col < BCS1.degree: is_pivot_1 = 0 is_pivot_2 = 0 - for i from cur_row <= i < BCS.dimension: + for i from cur_row <= i < BCS1.dimension: if bitset_check(&basis_1[i], gamma_1[cur_col]): is_pivot_1 = 1 piv_loc_1 = i break - for i from cur_row <= i < BCS.dimension: + for i from cur_row <= i < BCS1.dimension: if bitset_check(&basis_2[i], gamma_2[cur_col]): is_pivot_2 = 1 piv_loc_2 = i @@ -724,7 +830,7 @@ bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row]) if bitset_check(&basis_2[i], gamma_2[cur_col]): bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row]) - for i from cur_row < i < BCS.dimension: + for i from cur_row < i < BCS1.dimension: if bitset_check(&basis_1[i], gamma_1[cur_col]): bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row]) if bitset_check(&basis_2[i], gamma_2[cur_col]): @@ -736,34 +842,35 @@ return bitset_check(&basis_2[i], gamma_2[cur_col]) - bitset_check(&basis_1[i], gamma_1[cur_col]) return 0 -cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S): +cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, object S1, object S2): """ - Compare gamma_1(B) and gamma_2(B). + Compare gamma_1(S1) and gamma_2(S2). - Return return -1 if gamma_1(B) < gamma_2(B), 0 if gamma_1(B) == gamma_2(B), - 1 if gamma_1(B) > gamma_2(B). (Just like the python \code{cmp}) function. + Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2), + 1 if gamma_1(S1) > gamma_2(S2). (Just like the python \code{cmp}) function. INPUT: gamma_1, gamma_2 -- list permutations (inverse) - S -- a binary code struct object + S1, S2 -- a binary code struct object """ cdef int side=0, i, start, end, n_one_1, n_one_2, cur_col cdef int where_0, where_1 - cdef NonlinearBinaryCodeStruct BCS = S - cdef bitset_s *B_1_0 = BCS.scratch_bitsets # nwords of len degree - cdef bitset_s *B_1_1 = &BCS.scratch_bitsets[BCS.nwords] # nwords of len degree - cdef bitset_s *B_2_0 = &BCS.scratch_bitsets[2*BCS.nwords] # nwords of len degree - cdef bitset_s *B_2_1 = &BCS.scratch_bitsets[3*BCS.nwords] # nwords of len degree - cdef bitset_s *dividers = &BCS.scratch_bitsets[4*BCS.nwords] # 1 of len nwords + cdef NonlinearBinaryCodeStruct BCS1 = S1 + cdef NonlinearBinaryCodeStruct BCS2 = S2 + cdef bitset_s *B_1_0 = BCS1.scratch_bitsets # nwords of len degree + cdef bitset_s *B_1_1 = &BCS1.scratch_bitsets[BCS1.nwords] # nwords of len degree + cdef bitset_s *B_2_0 = &BCS1.scratch_bitsets[2*BCS1.nwords] # nwords of len degree + cdef bitset_s *B_2_1 = &BCS1.scratch_bitsets[3*BCS1.nwords] # nwords of len degree + cdef bitset_s *dividers = &BCS1.scratch_bitsets[4*BCS1.nwords] # 1 of len nwords cdef bitset_s *B_1_this, *B_1_other, *B_2_this, *B_2_other - for i from 0 <= i < BCS.nwords: - bitset_copy(&B_1_0[i], &BCS.words[i]) - bitset_copy(&B_2_0[i], &BCS.words[i]) + for i from 0 <= i < BCS1.nwords: + bitset_copy(&B_1_0[i], &BCS1.words[i]) + bitset_copy(&B_2_0[i], &BCS2.words[i]) bitset_zero(dividers) - bitset_set(dividers, BCS.nwords-1) + bitset_set(dividers, BCS1.nwords-1) - for cur_col from 0 <= cur_col < BCS.degree: + for cur_col from 0 <= cur_col < BCS1.degree: if side == 0: B_1_this = B_1_0 B_1_other = B_1_1 @@ -776,7 +883,7 @@ B_2_other = B_2_0 side ^= 1 start = 0 - while start < BCS.nwords: + while start < BCS1.nwords: end = start while not bitset_check(dividers, end): end += 1 @@ -960,7 +1067,7 @@ def random_tests(t=10.0, n_max=50, k_max=6, nwords_max=200, perms_per_code=10, density_range=(.1,.9)): """ Tests to make sure that C(gamma(B)) == C(B) for random permutations gamma - and random codes B. + and random codes B, and that is_isomorphic returns an isomorphism. INPUT: t -- run tests for approximately this many seconds @@ -977,7 +1084,8 @@ For each code B (B2) generated, we uniformly generate perms_per_code random permutations and verify that the canonical labels of B and the image of B - under the generated permutation are equal. + under the generated permutation are equal, and check that the double coset + function returns an isomorphism. DOCTEST: sage: import sage.groups.perm_gps.partn_ref.refinement_binary @@ -992,7 +1100,7 @@ from sage.combinat.permutation import Permutations from sage.matrix.constructor import random_matrix, matrix from sage.rings.finite_field import FiniteField as GF - cdef int h, i, j, n, k, num_tests = 0, num_codes = 0, passed = 1 + cdef int h, i, j, n, k, num_tests = 0, num_codes = 0 cdef LinearBinaryCodeStruct B, C cdef NonlinearBinaryCodeStruct B_n, C_n t_0 = walltime() @@ -1040,21 +1148,37 @@ B_n_M[j,B_n_relab[h]] = bitset_check(&B_n.words[j], h) C_n_M[j,C_n_relab[h]] = bitset_check(&C_n.words[j], h) if B_M.row_space() != C_M.row_space(): - print "B:" + print "can_lab error -- B:" for j from 0 <= j < B.dimension: print bitset_string(&B.basis[j]) print perm return if sorted(B_n_M.rows()) != sorted(C_n_M.rows()): - print "B_n:" + print "can_lab error -- B_n:" for j from 0 <= j < B_n.nwords: print bitset_string(&B_n.words[j]) print perm return + isom = B.is_isomorphic(C) + if not isom: + print "isom -- B:" + for j from 0 <= j < B.dimension: + print bitset_string(&B.basis[j]) + print perm + print isom + return + isom = B_n.is_isomorphic(C_n) + if not isom: + print "isom -- B_n:" + for j from 0 <= j < B_n.nwords: + print bitset_string(&B_n.words[j]) + print perm + print isom + return - num_tests += perms_per_code + num_tests += 4*perms_per_code num_codes += 2 - if passed: - print "All passed: %d random tests on %d codes."%(num_tests, num_codes) + + print "All passed: %d random tests on %d codes."%(num_tests, num_codes) diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_graphs.pxd --- a/sage/groups/perm_gps/partn_ref/refinement_graphs.pxd Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/refinement_graphs.pxd Mon Sep 15 17:34:52 2008 -0700 @@ -8,12 +8,14 @@ include '../../../ext/cdefs.pxi' include '../../../ext/stdsage.pxi' +include 'data_structures_pxd.pxi' # includes bitsets from sage.graphs.base.c_graph cimport CGraph from sage.graphs.base.sparse_graph cimport SparseGraph from sage.graphs.base.dense_graph cimport DenseGraph from sage.rings.integer cimport Integer -from sage.groups.perm_gps.partn_ref.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 +from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return +from double_coset cimport double_coset cdef class GraphStruct: cdef CGraph G @@ -22,7 +24,7 @@ cdef int *scratch # length 3n+1 cdef int refine_by_degree(PartitionStack *, object, int *, int) -cdef int compare_graphs(int *, int *, object) +cdef int compare_graphs(int *, int *, object, object) cdef bint all_children_are_equivalent(PartitionStack *, object) cdef inline int degree(PartitionStack *, CGraph, int, int, bint) cdef inline int sort_by_function(PartitionStack *, int, int *) diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_graphs.pyx --- a/sage/groups/perm_gps/partn_ref/refinement_graphs.pyx Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/refinement_graphs.pyx Mon Sep 15 17:34:52 2008 -0700 @@ -17,6 +17,144 @@ # Distributed under the terms of the GNU General Public License (GPL) # http://www.gnu.org/licenses/ #***************************************************************************** + +include 'data_structures_pyx.pxi' # includes bitsets + +def isomorphic(G1, G2, partition, ordering2, dig, use_indicator_function, sparse=False): + """ + Tests whether two graphs are isomorphic. + + sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import isomorphic + + sage: G = Graph(2) + sage: H = Graph(2) + sage: isomorphic(G, H, [[0,1]], [0,1], 0, 1) + [0, 1] + sage: isomorphic(G, H, [[0,1]], [0,1], 0, 1) + [0, 1] + sage: isomorphic(G, H, [[0],[1]], [0,1], 0, 1) + [0, 1] + sage: isomorphic(G, H, [[0],[1]], [1,0], 0, 1) + [1, 0] + + sage: G = Graph(3) + sage: H = Graph(3) + sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1) + [0, 1, 2] + sage: G.add_edge(0,1) + sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1) + False + sage: H.add_edge(1,2) + sage: isomorphic(G, H, [[0,1,2]], [0,1,2], 0, 1) + [1, 2, 0] + + """ + cdef int **part + cdef int *output, *ordering + cdef CGraph G + cdef GraphStruct GS1 = GraphStruct() + cdef GraphStruct GS2 = GraphStruct() + cdef GraphStruct GS + cdef int i, j, n = -1 + + from sage.graphs.graph import GenericGraph, Graph, DiGraph + for G_in in [G1, G2]: + if G_in is G1: + GS = GS1 + else: + GS = GS2 + if isinstance(G_in, GenericGraph): + if n == -1: + n = G_in.num_verts() + elif n != G_in.num_verts(): + return False + G_in = G_in.copy() + if G_in.vertices() != range(n): + to = G_in.relabel(return_map=True) + frm = {} + for v in to.iterkeys(): + frm[to[v]] = v + partition = [[to[v] for v in cell] for cell in partition] + else: + to = range(n) + frm = to + if sparse: + G = SparseGraph(n) + else: + G = DenseGraph(n) + if G_in.is_directed(): + for i from 0 <= i < n: + for _,j,_ in G_in.outgoing_edge_iterator(i): + G.add_arc(i,j) + else: + for i from 0 <= i < n: + for _,j,_ in G_in.edge_iterator(i): + if j <= i: + G.add_arc(i,j) + G.add_arc(j,i) + elif isinstance(G_in, CGraph): + G = G_in + if n == -1: + n = G.num_verts + elif n != G.num_verts: + return False + to = range(n) + frm = to + else: + raise TypeError("G must be a Sage graph.") + GS.G = G + GS.directed = 1 if dig else 0 + GS.use_indicator = 1 if use_indicator_function else 0 + + if n == 0: + return {} + + part = sage_malloc((len(partition)+1) * sizeof(int *)) + ordering = sage_malloc(n * sizeof(int)) + if part is NULL or ordering is NULL: + if part is not NULL: sage_free(part) + if ordering is not NULL: sage_free(ordering) + raise MemoryError + for i from 0 <= i < len(partition): + part[i] = sage_malloc((len(partition[i])+1) * sizeof(int)) + if part[i] is NULL: + for j from 0 <= j < i: + sage_free(part[j]) + sage_free(part) + raise MemoryError + for j from 0 <= j < len(partition[i]): + part[i][j] = partition[i][j] + part[i][len(partition[i])] = -1 + part[len(partition)] = NULL + for i from 0 <= i < n: + ordering[i] = ordering2[i] + + GS1.scratch = sage_malloc((3*n+1) * sizeof(int)) + GS2.scratch = sage_malloc((3*n+1) * sizeof(int)) + if GS1.scratch is NULL or GS2.scratch is NULL: + if GS1.scratch is not NULL: sage_free(GS1.scratch) + if GS2.scratch is not NULL: sage_free(GS2.scratch) + for j from 0 <= j < len(partition): + sage_free(part[j]) + sage_free(part) + raise MemoryError + + output = double_coset(GS1, GS2, part, ordering, n, &all_children_are_equivalent, &refine_by_degree, &compare_graphs) + + for i from 0 <= i < len(partition): + sage_free(part[i]) + sage_free(part) + sage_free(ordering) + sage_free(GS1.scratch) + sage_free(GS2.scratch) + + if output is NULL: + return False + else: + output_py = [output[i] for i from 0 <= i < n] + sage_free(output) + # TODO: figure out frm, to stuff to relabel for consistency with input + return output_py def search_tree(G_in, partition, lab=True, dig=False, dict_rep=False, certify=False, verbosity=0, use_indicator_function=True, sparse=False, @@ -200,6 +338,19 @@ sage: st(G, [range(G.num_verts())])[1] == st(H, [range(H.num_verts())])[1] True + sage: from sage.graphs.graph import graph_isom_equivalent_non_multi_graph + sage: G = Graph(multiedges=True, implementation='networkx') + sage: G.add_edge(('a', 'b')) + sage: G.add_edge(('a', 'b')) + sage: G.add_edge(('a', 'b')) + sage: G, Pi = graph_isom_equivalent_non_multi_graph(G, [['a','b']]) + sage: s,b = st(G, Pi, lab=False, dict_rep=True) + sage: sorted(b.items()) + [(('o', 'a'), 5), (('o', 'b'), 1), (('x', 0), 2), (('x', 1), 3), (('x', 2), 4)] + + sage: st(Graph(':Dkw'), [range(5)], lab=False, dig=True) + [[4, 1, 2, 3, 0], [0, 2, 1, 3, 4]] + """ cdef CGraph G cdef int i, j, n @@ -277,10 +428,7 @@ if dict_rep: ddd = {} for v in frm.iterkeys(): - if frm[v] != 0: - ddd[v] = frm[v] - else: - ddd[v] = n + ddd[frm[v]] = v if v != 0 else n return_tuple.append(ddd) if lab: if isinstance(G_in, GenericGraph): @@ -348,6 +496,7 @@ cdef int current_cell, i, r cdef int first_largest_subcell cdef int invariant = 1 + cdef int max_degree cdef int *degrees = GS.scratch # length 3n+1 cdef bint necessary_to_split_cell cdef int against_index @@ -358,10 +507,13 @@ invariant += 50 i = current_cell necessary_to_split_cell = 0 + max_degree = 0 while 1: degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 0) if degrees[i-current_cell] != degrees[0]: necessary_to_split_cell = 1 + if degrees[i-current_cell] > max_degree: + max_degree = degrees[i-current_cell] i += 1 if PS.levels[i-1] <= PS.depth: break @@ -369,7 +521,7 @@ if necessary_to_split_cell: invariant += 10 first_largest_subcell = sort_by_function(PS, current_cell, degrees) - invariant += first_largest_subcell + invariant += first_largest_subcell + max_degree against_index = current_cell_against while against_index < ctrb_len: if cells_to_refine_by[against_index] == current_cell: @@ -385,7 +537,6 @@ r += 1 if r >= i: break - invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 0) invariant += (i - current_cell) current_cell = i if GS.directed: @@ -396,10 +547,13 @@ invariant += 20 i = current_cell necessary_to_split_cell = 0 + max_degree = 0 while 1: degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 1) if degrees[i-current_cell] != degrees[0]: necessary_to_split_cell = 1 + if degrees[i-current_cell] > max_degree: + max_degree = degrees[i-current_cell] i += 1 if PS.levels[i-1] <= PS.depth: break @@ -407,7 +561,7 @@ if necessary_to_split_cell: invariant += 7 first_largest_subcell = sort_by_function(PS, current_cell, degrees) - invariant += first_largest_subcell + invariant += first_largest_subcell + max_degree against_index = current_cell_against while against_index < ctrb_len: if cells_to_refine_by[against_index] == current_cell: @@ -425,7 +579,6 @@ r += 1 if r >= i: break - invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 1) invariant += (i - current_cell) current_cell = i current_cell_against += 1 @@ -434,28 +587,30 @@ else: return 0 -cdef int compare_graphs(int *gamma_1, int *gamma_2, object S): +cdef int compare_graphs(int *gamma_1, int *gamma_2, object S1, object S2): r""" - Compare gamma_1(G) and gamma_2(G). + Compare gamma_1(S1) and gamma_2(S2). - Return return -1 if gamma_1(G) < gamma_2(G), 0 if gamma_1(G) == - gamma_2(G), 1 if gamma_1(G) > gamma_2(G). (Just like the python + Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == + gamma_2(S2), 1 if gamma_1(S1) > gamma_2(S2). (Just like the python \code{cmp}) function. INPUT: gamma_1, gamma_2 -- list permutations (inverse) - S -- a graph struct object + S1, S2 -- graph struct objects """ cdef int i, j - cdef GraphStruct GS = S - cdef CGraph G = GS.G - for i from 0 <= i < G.num_verts: - for j from 0 <= j < G.num_verts: - if G.has_arc_unsafe(gamma_1[i], gamma_1[j]): - if not G.has_arc_unsafe(gamma_2[i], gamma_2[j]): + cdef GraphStruct GS1 = S1 + cdef GraphStruct GS2 = S2 + cdef CGraph G1 = GS1.G + cdef CGraph G2 = GS2.G + for i from 0 <= i < G1.num_verts: + for j from 0 <= j < G1.num_verts: + if G1.has_arc_unsafe(gamma_1[i], gamma_1[j]): + if not G2.has_arc_unsafe(gamma_2[i], gamma_2[j]): return 1 - elif G.has_arc_unsafe(gamma_2[i], gamma_2[j]): + elif G2.has_arc_unsafe(gamma_2[i], gamma_2[j]): return -1 return 0 @@ -472,6 +627,8 @@ S -- a graph struct object """ cdef GraphStruct GS = S + if GS.directed: + return 0 cdef CGraph G = GS.G cdef int i, n = PS.degree cdef bint in_cell = 0 @@ -622,7 +779,7 @@ def random_tests(t=10.0, n_max=60, perms_per_graph=10): """ Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma - and random graphs G. + and random graphs G, and that isomorphic returns an isomorphism. INPUT: t -- run tests for approximately this many seconds @@ -637,7 +794,8 @@ For each graph G generated, we uniformly generate perms_per_graph random permutations and verify that the canonical labels of G and the image of G - under the generated permutation are equal. + under the generated permutation are equal, and that the isomorphic function + returns an isomorphism. TESTS: sage: import sage.groups.perm_gps.partn_ref.refinement_graphs @@ -651,7 +809,7 @@ from sage.misc.prandom import random, randint from sage.graphs.graph_generators import GraphGenerators, DiGraphGenerators from sage.combinat.permutation import Permutations - cdef int i, j, num_tests = 0, num_graphs = 0, passed = 1 + cdef int i, j, num_tests = 0, num_graphs = 0 GG = GraphGenerators() DGG = DiGraphGenerators() t_0 = walltime() @@ -663,16 +821,23 @@ G = GG.RandomGNP(n, p) H = G.copy() for i from 0 <= i < perms_per_graph: + G = H.copy() G1 = search_tree(G, [G.vertices()])[1] perm = list(S.random_element()) perm = [perm[j]-1 for j from 0 <= j < n] G.relabel(perm) G2 = search_tree(G, [G.vertices()])[1] if G1 != G2: - print "FAILURE: graph6-" + print "search_tree FAILURE: graph6-" print H.graph6_string() print perm - passed = 0 + return + isom = isomorphic(G, H, [range(n)], range(n), 0, 1) + if not isom or G.relabel(isom, inplace=False) != H: + print "isom FAILURE: graph6-" + print H.graph6_string() + print perm + return D = DGG.RandomDirectedGNP(n, p) D.loops(True) @@ -681,19 +846,26 @@ D.add_edge(i,i) E = D.copy() for i from 0 <= i < perms_per_graph: + D = E.copy() D1 = search_tree(D, [D.vertices()], dig=True)[1] perm = list(S.random_element()) perm = [perm[j]-1 for j from 0 <= j < n] D.relabel(perm) D2 = search_tree(D, [D.vertices()], dig=True)[1] if D1 != D2: - print "FAILURE: dig6-" + print "search_tree FAILURE: dig6-" print E.dig6_string() print perm - passed = 0 - num_tests += 2*perms_per_graph + return + isom = isomorphic(D, E, [range(n)], range(n), 1, 1) + if not isom or D.relabel(isom, inplace=False) != E: + print "isom FAILURE: dig6-" + print E.dig6_string() + print perm + print isom + return + num_tests += 4*perms_per_graph num_graphs += 2 - if passed: - print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs) + print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs) diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_matrices.pxd --- a/sage/groups/perm_gps/partn_ref/refinement_matrices.pxd Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/refinement_matrices.pxd Mon Sep 15 17:34:52 2008 -0700 @@ -8,12 +8,13 @@ include '../../../ext/cdefs.pxi' include '../../../ext/stdsage.pxi' -include '../../../misc/bitset_pxd.pxi' +include 'data_structures_pxd.pxi' # includes bitsets from sage.rings.integer cimport Integer -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 -from sage.groups.perm_gps.partn_ref.refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree -from sage.groups.perm_gps.partn_ref.refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent +from automorphism_group_canonical_label cimport get_aut_gp_and_can_lab, aut_gp_and_can_lab_return +from refinement_binary cimport NonlinearBinaryCodeStruct, refine_by_bip_degree +from refinement_binary cimport all_children_are_equivalent as all_binary_children_are_equivalent +from double_coset cimport double_coset cdef class MatrixStruct: cdef list symbol_structs @@ -26,6 +27,6 @@ cdef aut_gp_and_can_lab_return *output cdef int refine_matrix(PartitionStack *, object, int *, int) -cdef int compare_matrices(int *, int *, object) +cdef int compare_matrices(int *, int *, object, object) cdef bint all_matrix_children_are_equivalent(PartitionStack *, object) diff -r 0d158c58d1c2 -r 2df5cad14396 sage/groups/perm_gps/partn_ref/refinement_matrices.pyx --- a/sage/groups/perm_gps/partn_ref/refinement_matrices.pyx Mon Sep 15 05:47:46 2008 -0700 +++ b/sage/groups/perm_gps/partn_ref/refinement_matrices.pyx Mon Sep 15 17:34:52 2008 -0700 @@ -22,7 +22,8 @@ # http://www.gnu.org/licenses/ #***************************************************************************** -include '../../../misc/bitset.pxi' +include 'data_structures_pyx.pxi' # includes bitsets + from sage.misc.misc import uniq from sage.matrix.constructor import Matrix @@ -127,6 +128,11 @@ """ cdef int **part, i, j + cdef NonlinearBinaryCodeStruct S_temp + for i from 0 <= i < self.nsymbols: + S_temp = self.symbol_structs[i] + S_temp.first_time = 1 + if partition is None: partition = [range(self.degree)] part = sage_malloc((len(partition)+1) * sizeof(int *)) @@ -195,6 +201,62 @@ self.run() return [self.output.relabeling[i] for i from 0 <= i < self.degree] + def is_isomorphic(self, MatrixStruct other): + """ + Calculate whether self is isomorphic to other. + + EXAMPLES: + sage: from sage.groups.perm_gps.partn_ref.refinement_matrices import MatrixStruct + sage: M = MatrixStruct(Matrix(GF(11), [[1,2,3,0,0,0],[0,0,0,1,2,3]])) + sage: N = MatrixStruct(Matrix(GF(11), [[0,1,0,2,0,3],[1,0,2,0,3,0]])) + sage: M.is_isomorphic(N) + [0, 2, 4, 1, 3, 5] + + """ + + cdef int **part, i, j + cdef int *output, *ordering + cdef NonlinearBinaryCodeStruct S_temp + for i from 0 <= i < self.nsymbols: + S_temp = self.symbol_structs[i] + S_temp.first_time = 1 + S_temp = other.symbol_structs[i] + S_temp.first_time = 1 + partition = [range(self.degree)] + part = sage_malloc((len(partition)+1) * sizeof(int *)) + ordering = sage_malloc(self.degree * sizeof(int)) + if part is NULL or ordering is NULL: + if part is not NULL: sage_free(part) + if ordering is not NULL: sage_free(ordering) + raise MemoryError + for i from 0 <= i < len(partition): + part[i] = sage_malloc((len(partition[i])+1) * sizeof(int)) + if part[i] is NULL: + for j from 0 <= j < i: + sage_free(part[j]) + sage_free(part) + raise MemoryError + for j from 0 <= j < len(partition[i]): + part[i][j] = partition[i][j] + part[i][len(partition[i])] = -1 + part[len(partition)] = NULL + for i from 0 <= i < self.degree: + ordering[i] = i + + output = double_coset(self, other, part, ordering, self.degree, &all_matrix_children_are_equivalent, &refine_matrix, &compare_matrices) + + for i from 0 <= i < len(partition): + sage_free(part[i]) + sage_free(part) + sage_free(ordering) + + if output is NULL: + return False + else: + output_py = [output[i] for i from 0 <= i < self.degree] + sage_free(output) + return output_py + cdef int refine_matrix(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len): cdef MatrixStruct M = S cdef int i, temp_inv, invariant = 1 @@ -209,16 +271,18 @@ changed = 0 return invariant -cdef int compare_matrices(int *gamma_1, int *gamma_2, object S): - cdef MatrixStruct MS = S - M = MS.matrix +cdef int compare_matrices(int *gamma_1, int *gamma_2, object S1, object S2): + cdef MatrixStruct MS1 = S1 + cdef MatrixStruct MS2 = S2 + M1 = MS1.matrix + M2 = MS2.matrix cdef int i - M1 = Matrix(M.base_ring(), M.nrows(), M.ncols(), sparse=M.is_sparse()) - M2 = Matrix(M.base_ring(), M.nrows(), M.ncols(), sparse=M.is_sparse()) - for i from 0 <= i < M.ncols(): - M1.set_column(i, M.column(gamma_1[i])) - M2.set_column(i, M.column(gamma_2[i])) - return cmp(sorted(M1.rows()), sorted(M2.rows())) + MM1 = Matrix(M1.base_ring(), M1.nrows(), M1.ncols(), sparse=M1.is_sparse()) + MM2 = Matrix(M2.base_ring(), M2.nrows(), M2.ncols(), sparse=M2.is_sparse()) + for i from 0 <= i < M1.ncols(): + MM1.set_column(i, M1.column(gamma_1[i])) + MM2.set_column(i, M2.column(gamma_2[i])) + return cmp(sorted(MM1.rows()), sorted(MM2.rows())) cdef bint all_matrix_children_are_equivalent(PartitionStack *PS, object S): cdef MatrixStruct M = S @@ -230,7 +294,8 @@ def random_tests(t=10.0, nrows_max=50, ncols_max=50, nsymbols_max=20, perms_per_matrix=10, density_range=(.1,.9)): """ Tests to make sure that C(gamma(M)) == C(M) for random permutations gamma - and random matrices M. + and random matrices M, and that M.is_isomorphic(gamma(M)) returns an + isomorphism. INPUT: t -- run tests for approximately this many seconds @@ -247,7 +312,8 @@ For each matrix M generated, we uniformly generate perms_per_matrix random permutations and verify that the canonical labels of M and the image of M - under the generated permutation are equal. + under the generated permutation are equal, and that the isomorphism is + discovered by the double coset function. DOCTEST: sage: import sage.groups.perm_gps.partn_ref.refinement_matrices @@ -263,7 +329,7 @@ from sage.matrix.constructor import random_matrix, matrix from sage.rings.finite_field import FiniteField as GF from sage.rings.arith import next_prime - cdef int h, i, j, nrows, k, num_tests = 0, num_matrices = 0, passed = 1 + cdef int h, i, j, nrows, k, num_tests = 0, num_matrices = 0 cdef MatrixStruct M, N t_0 = walltime() while walltime(t_0) < t: @@ -300,15 +366,22 @@ if M_C != N_C: print "M:" - print M + print M.matrix.str() print "perm:" print perm + return + + isom = M.is_isomorphic(N) + if not isom: + print "isom FAILURE: M:" + print M.matrix.str() + print "isom FAILURE: N:" + print N.matrix.str() return num_tests += perms_per_matrix num_matrices += 2 - if passed: - print "All passed: %d random tests on %d matrices."%(num_tests, num_matrices) + print "All passed: %d random tests on %d matrices."%(num_tests, num_matrices) diff -r 0d158c58d1c2 -r 2df5cad14396 setup.py --- a/setup.py Mon Sep 15 05:47:46 2008 -0700 +++ b/setup.py Mon Sep 15 17:34:52 2008 -0700 @@ -717,6 +717,9 @@ Extension('sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label', sources = ['sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx']), \ + + Extension('sage.groups.perm_gps.partn_ref.double_coset', + sources = ['sage/groups/perm_gps/partn_ref/double_coset.pyx']), \ Extension('sage.groups.perm_gps.partn_ref.refinement_graphs', sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx']), \