# HG changeset patch # User Robert L. Miller # Date 1216425781 25200 # Node ID 074731a2698613667ad1b649a96a8446a892ffd6 # Parent 91af4c3f6b92316e4117af342d82be33cf5d949b Refactored graph isomorphism code Split up bitset.pxi to allow use in pxd files. diff -r 5e2426d277eb sage/graphs/graph_generators.py --- a/sage/graphs/graph_generators.py Wed Jul 30 06:34:58 2008 -0700 +++ b/sage/graphs/graph_generators.py Mon Aug 11 18:26:07 2008 -0700 @@ -2991,6 +2991,35 @@ class DiGraphGenerators(): import networkx return graph.DiGraph(networkx.gnc_graph(n, seed)) + def RandomDirectedGNP(self, n, p): + r""" + Returns a random digraph on $n$ nodes. Each edge is inserted + independently with probability $p$. + + REFERENCES: + [1] P. Erdos and A. Renyi, On Random Graphs, Publ. Math. 6, 290 (1959). + [2] E. N. Gilbert, Random Graphs, Ann. Math. Stat., 30, 1141 (1959). + + PLOTTING: + When plotting, this graph will use the default spring-layout + algorithm, unless a position dictionary is specified. + + EXAMPLE: + sage: digraphs.RandomDirectedGNP(10, .2).num_verts() + 10 + + """ + from random import random + D = graph.DiGraph(n) + for i in xrange(n): + for j in xrange(i): + if random() < p: + D.add_edge(i,j) + for j in xrange(i+1,n): + if random() < p: + D.add_edge(i,j) + return D + def RandomDirectedGNR(self, n, p, seed=None): """ Returns a random GNR (growing network with redirection) digraph with n diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/__init__.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/__init__.py Mon Aug 11 18:26:07 2008 -0700 @@ -0,0 +1,1 @@ + diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd Mon Aug 11 18:26:07 2008 -0700 @@ -0,0 +1,75 @@ + +#***************************************************************************** +# 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' + +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 + int num_gens + int *relabeling + int *base + int base_size + mpz_t order + +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) + + + diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx Mon Aug 11 18:26:07 2008 -0700 @@ -0,0 +1,1324 @@ +r""" +Partition refinement algorithms for computing automorphism groups and canonical +labels of various objects. + +This module implements a general algorithm for computing automorphism groups and +canonical labels of objects. The class of objects in question must be some kind +of structure for which an automorphism is a permutation in $S_n$ for some $n$, +which we call here the order of the object. It should be noted that the word +"canonical" in the term canonical label is used loosely. Given an object $X$, +the canonical label $C(X)$ is just an object for which $C(X) = C(Y)$ for any +$Y$ such that $X \cong Y$. + +To understand the algorithm better, a few definitions will help. First one +should note that throughout this module, $S_n$ is defined as the set of +permutations of the set $[n] := \{0, 1, ..., n-1\}$. One speaks of +ordered partitions $\Pi = (C_1, ..., C_k)$ of $[n]$. The $C_i$ are disjoint +subsets of $[n]$ whose union is equal to $[n]$, and we call them cells. A +partition $\Pi_1$ is said to be finer than another partition $\Pi_2$ if every +cell of $\Pi_1$ is a subset of some cell of $\Pi_2$ (in this situation we also +say that $\Pi_2$ is coarser than $\Pi_1$). A partition stack +$\nu = (\Pi_0, ..., \Pi_k)$ is a sequence of partitions such that $\Pi_i$ is +finer than $\Pi_{i-1}$ for each $i$ such that $1 \leq i \leq k$. Sometimes these +are called nested partitions. The depth of $\nu$ is defined to be $k$ and the +degree of $\nu$ is defined to be $n$. Partition stacks are implemented as +\code{PartitionStack} (see automorphism_group_canonical_label.pxd for more +information). Finally, we say that a permutation respects the partition $\Pi$ if +its orbits form a partition which is finer than $\Pi$. + +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{get_aut_gp_and_can_lab}. 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.automorphism_group_canonical_label + +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 + +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 + +cdef int refine_and_return_invariant_trivial(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len): + return 0 + +cdef int compare_structures_trivial(int *gamma_1, int *gamma_2, object S): + return 0 + +def test_get_aut_gp_and_can_lab_trivially(int n=6, partition=[[0,1,2],[3,4],[5]], + canonical_label=True, base=False): + """ + sage: tttt = sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label.test_get_aut_gp_and_can_lab_trivially + sage: tttt() + 12 + sage: tttt(canonical_label=False, base=False) + 12 + sage: tttt(canonical_label=False, base=True) + 12 + sage: tttt(canonical_label=True, base=True) + 12 + sage: tttt(n=0, partition=[]) + 1 + sage: tttt(n=0, partition=[], canonical_label=False, base=False) + 1 + sage: tttt(n=0, partition=[], canonical_label=False, base=True) + 1 + sage: tttt(n=0, partition=[], canonical_label=True, base=True) + 1 + + """ + cdef int i, j + cdef aut_gp_and_can_lab_return *output + cdef Integer I = Integer(0) + cdef int **part = sage_malloc((len(partition)+1) * sizeof(int *)) + for i from 0 <= i < len(partition): + part[i] = sage_malloc((len(partition[i])+1) * sizeof(int)) + for j from 0 <= j < len(partition[i]): + part[i][j] = partition[i][j] + part[i][len(partition[i])] = -1 + part[len(partition)] = NULL + output = get_aut_gp_and_can_lab([], part, n, &all_children_are_equivalent_trivial, &refine_and_return_invariant_trivial, &compare_structures_trivial, canonical_label, base, True) + mpz_set(I.value, output.order) + print I + for i from 0 <= i < len(partition): + sage_free(part[i]) + sage_free(part) + mpz_clear(output.order) + sage_free(output.generators) + if base: + sage_free(output.base) + if canonical_label: + sage_free(output.relabeling) + sage_free(output) + +cdef aut_gp_and_can_lab_return *get_aut_gp_and_can_lab(object S, int **partition, int n, + bint (*all_children_are_equivalent)(PartitionStack *PS, object S), + int (*refine_and_return_invariant)\ + (PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len), + int (*compare_structures)(int *gamma_1, int *gamma_2, object S), + bint canonical_label, bint base, bint order): + """ + Traverse the search space for subgroup/canonical label calculation. + + INPUT: + S -- pointer to the structure + partition -- array representing a partition of the points + len_partition -- length of the partition + n -- the number of points (points are assumed to be 0,1,...,n-1) + canonical_label -- whether to search for canonical label - if True, return + the permutation taking S to its canonical label + base -- if True, return the base for which the given generating set is a + strong generating set + order -- if True, return the order of the automorphism group + all_children_are_equivalent -- pointer to a function + INPUT: + PS -- pointer to a partition stack + S -- pointer to the structure + OUTPUT: + bint -- returns True if it can be determined that all refinements below + the current one will result in an equivalent discrete partition + refine_and_return_invariant -- pointer to a function + INPUT: + PS -- pointer to a partition stack + S -- pointer to the structure + alpha -- an array consisting of numbers, which indicate the starting + positions of the cells to refine against (will likely be modified) + OUTPUT: + int -- returns an invariant under application of arbitrary permutations + compare_structures -- pointer to a function + INPUT: + gamma_1, gamma_2 -- (list) permutations of the points of S + S -- pointer to the structure + OUTPUT: + int -- 0 if gamma_1(S) = gamma_2(S), otherwise -1 or 1 (see docs for cmp), + such that the set of all structures is well-ordered + OUTPUT: + pointer to a aut_gp_and_can_lab_return struct + + """ + cdef PartitionStack *current_ps, *first_ps, *label_ps + cdef int first_meets_current = -1 + cdef int label_meets_current + cdef int current_kids_are_same = 1 + cdef int first_kids_are_same + + cdef int *current_indicators, *first_indicators, *label_indicators + cdef int first_and_current_indicator_same + cdef int label_and_current_indicator_same = -1 + cdef int compared_current_and_label_indicators + + cdef OrbitPartition *orbits_of_subgroup + cdef mpz_t subgroup_size + cdef int subgroup_primary_orbit_size = 0 + cdef int minimal_in_primary_orbit + cdef int size_of_generator_array = 2*n*n + + cdef bitset_t *fixed_points_of_generators # i.e. fp + cdef bitset_t *minimal_cell_reps_of_generators # i.e. mcr + cdef int len_of_fp_and_mcr = 100 + cdef int index_in_fp_and_mcr = -1 + + cdef bitset_t *vertices_to_split + cdef bitset_t vertices_have_been_reduced + cdef int *permutation, *id_perm, *cells_to_refine_by + cdef int *vertices_determining_current_stack + + cdef int i, j, k + cdef bint discrete, automorphism, update_label + cdef bint backtrack, new_vertex, narrow, mem_err = 0 + + cdef aut_gp_and_can_lab_return *output = sage_malloc(sizeof(aut_gp_and_can_lab_return)) + if output is NULL: + raise MemoryError + + if n == 0: + output.generators = NULL + output.num_gens = 0 + output.relabeling = NULL + output.base = NULL + output.base_size = 0 + mpz_init_set_si(output.order, 1) + return output + + # Allocate: + mpz_init_set_si(subgroup_size, 1) + + output.generators = sage_malloc( size_of_generator_array * sizeof(int) ) + output.num_gens = 0 + if base: + output.base = sage_malloc( n * sizeof(int) ) + output.base_size = 0 + + current_indicators = sage_malloc( n * sizeof(int) ) + first_indicators = sage_malloc( n * sizeof(int) ) + + if canonical_label: + label_indicators = sage_malloc( n * sizeof(int) ) + output.relabeling = 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) ) + + current_ps = PS_new(n, 0) + orbits_of_subgroup = OP_new(n) + + # Check for allocation failures: + cdef bint succeeded = 1 + if canonical_label: + if label_indicators is NULL or output.relabeling is NULL: + succeeded = 0 + if label_indicators is not NULL: + sage_free(label_indicators) + if output.relabeling is not NULL: + sage_free(output.relabeling) + if base: + if output.base is NULL: + succeeded = 0 + if not succeeded or \ + output.generators is NULL or \ + current_indicators is NULL or \ + first_indicators is NULL or \ + fixed_points_of_generators is NULL or \ + minimal_cell_reps_of_generators is NULL or \ + vertices_to_split is NULL or \ + permutation is NULL or \ + id_perm is NULL or \ + cells_to_refine_by is NULL or \ + vertices_determining_current_stack is NULL or \ + current_ps is NULL or \ + orbits_of_subgroup is NULL: + if output.generators is not NULL: + sage_free(output.generators) + if base: + if output.base is not NULL: + sage_free(output.base) + if current_indicators is not NULL: + sage_free(current_indicators) + if first_indicators is not NULL: + sage_free(first_indicators) + if fixed_points_of_generators is not NULL: + sage_free(fixed_points_of_generators) + if minimal_cell_reps_of_generators is not NULL: + sage_free(minimal_cell_reps_of_generators) + if vertices_to_split is not NULL: + sage_free(vertices_to_split) + if permutation is not NULL: + sage_free(permutation) + if id_perm is not NULL: + sage_free(id_perm) + if cells_to_refine_by is not NULL: + sage_free(cells_to_refine_by) + if vertices_determining_current_stack is not NULL: + sage_free(vertices_determining_current_stack) + if current_ps is not NULL: + PS_dealloc(current_ps) + if orbits_of_subgroup is not NULL: + OP_dealloc(orbits_of_subgroup) + sage_free(output) + mpz_clear(subgroup_size) + raise MemoryError + + # Initialize bitsets, checking for allocation failures: + succeeded = 1 + for i from 0 <= i < len_of_fp_and_mcr: + try: + bitset_init(fixed_points_of_generators[i], n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < i: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + break + try: + bitset_init(minimal_cell_reps_of_generators[i], n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < i: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + bitset_clear(fixed_points_of_generators[i]) + break + if succeeded: + for i from 0 <= i < n: + try: + bitset_init(vertices_to_split[i], n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < i: + bitset_clear(vertices_to_split[j]) + for j from 0 <= j < len_of_fp_and_mcr: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + break + if succeeded: + try: + bitset_init(vertices_have_been_reduced, n) + except MemoryError: + succeeded = 0 + for j from 0 <= j < n: + bitset_clear(vertices_to_split[j]) + for j from 0 <= j < len_of_fp_and_mcr: + bitset_clear(fixed_points_of_generators[j]) + bitset_clear(minimal_cell_reps_of_generators[j]) + if not succeeded: + sage_free(output.generators) + if base: + sage_free(output.base) + sage_free(current_indicators) + sage_free(first_indicators) + if canonical_label: + sage_free(output.relabeling) + sage_free(label_indicators) + sage_free(fixed_points_of_generators) + sage_free(minimal_cell_reps_of_generators) + sage_free(vertices_to_split) + sage_free(permutation) + sage_free(id_perm) + sage_free(cells_to_refine_by) + sage_free(vertices_determining_current_stack) + PS_dealloc(current_ps) + sage_free(output) + mpz_clear(subgroup_size) + raise MemoryError + + bitset_zero(vertices_have_been_reduced) + + # Copy data of partition to current_ps. + i = 0 + j = 0 + while partition[i] is not NULL: + k = 0 + while partition[i][k] != -1: + current_ps.entries[j+k] = partition[i][k] + current_ps.levels[j+k] = n + k += 1 + current_ps.levels[j+k-1] = 0 + PS_move_min_to_front(current_ps, j, j+k-1) + j += k + i += 1 + current_ps.levels[j-1] = -1 + current_ps.depth = 0 + current_ps.degree = n + + # default values of "infinity" + for i from 0 <= i < n: + first_indicators[i] = -1 + if canonical_label: + for i from 0 <= i < n: + label_indicators[i] = -1 + + # set up the identity permutation + for i from 0 <= i < n: + id_perm[i] = i + + # Our first refinement needs to check every cell of the partition, + # so cells_to_refine_by needs to be a list of pointers to each cell. + j = 1 + cells_to_refine_by[0] = 0 + for i from 0 < i < n: + if current_ps.levels[i-1] == 0: + cells_to_refine_by[j] = i + j += 1 + # Ignore the invariant this time, since we are + # creating the root of the search tree. + refine_and_return_invariant(current_ps, S, cells_to_refine_by, j) + + # Refine down to a discrete partition + while not PS_is_discrete(current_ps): + if not all_children_are_equivalent(current_ps, S): + current_kids_are_same = current_ps.depth + 1 + vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth]) + bitset_unset(vertices_have_been_reduced, current_ps.depth) + i = current_ps.depth + current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by) + first_indicators[i] = current_indicators[i] + if canonical_label: + label_indicators[i] = current_indicators[i] + if base: + output.base[i] = vertices_determining_current_stack[i] + output.base_size += 1 + + first_meets_current = current_ps.depth + first_kids_are_same = current_ps.depth + first_and_current_indicator_same = current_ps.depth + first_ps = PS_copy(current_ps) + if canonical_label: + label_meets_current = current_ps.depth + label_and_current_indicator_same = current_ps.depth + compared_current_and_label_indicators = 0 + label_ps = PS_copy(current_ps) + current_ps.depth -= 1 + + # Main loop: + while current_ps.depth != -1: + + # If necessary, update label_meets_current + if canonical_label and label_meets_current > current_ps.depth: + label_meets_current = current_ps.depth + + # I. Search for a new vertex to split, and update subgroup information + new_vertex = 0 + if current_ps.depth > first_meets_current: + # If we are not at a node of the first stack, reduce size of + # vertices_to_split by using the symmetries we already know. + if not bitset_check(vertices_have_been_reduced, current_ps.depth): + for i from 0 <= i <= index_in_fp_and_mcr: + j = 0 + while j < current_ps.depth and bitset_check(fixed_points_of_generators[i], vertices_determining_current_stack[j]): + j += 1 + # If each vertex split so far is fixed by generator i, + # then remove elements of vertices_to_split which are + # not minimal in their orbits under generator i. + if j == current_ps.depth: + for k from 0 <= k < n: + if bitset_check(vertices_to_split[current_ps.depth], k) and not bitset_check(minimal_cell_reps_of_generators[i], k): + bitset_flip(vertices_to_split[current_ps.depth], k) + bitset_flip(vertices_have_been_reduced, current_ps.depth) + # Look for a new point to split. + i = vertices_determining_current_stack[current_ps.depth] + 1 + while i < n and not bitset_check(vertices_to_split[current_ps.depth], i): + i += 1 + if i < n: + # There is a new point. + vertices_determining_current_stack[current_ps.depth] = i + new_vertex = 1 + else: + # No new point: backtrack. + current_ps.depth -= 1 + else: + # If we are at a node of the first stack, the above reduction + # will not help. Also, we must update information about + # primary orbits here. + if current_ps.depth < first_meets_current: + # If we are done searching under this part of the first + # stack, then first_meets_current is one higher, and we + # are looking at a new primary orbit (corresponding to a + # larger subgroup in the stabilizer chain). + first_meets_current = current_ps.depth + for i from 0 <= i < n: + if bitset_check(vertices_to_split[current_ps.depth], i): + minimal_in_primary_orbit = i + break + while 1: + i = vertices_determining_current_stack[current_ps.depth] + # This was the last point to be split here. + # If it is in the same orbit as minimal_in_primary_orbit, + # then count it as an element of the primary orbit. + if OP_find(orbits_of_subgroup, i) == OP_find(orbits_of_subgroup, minimal_in_primary_orbit): + subgroup_primary_orbit_size += 1 + # Look for a new point to split. + i += 1 + while i < n and not bitset_check(vertices_to_split[current_ps.depth], i): + i += 1 + if i < n: + # There is a new point. + vertices_determining_current_stack[current_ps.depth] = i + if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, i)] == i: + new_vertex = 1 + break + else: + # No new point: backtrack. + # Note that now, we are backtracking up the first stack. + vertices_determining_current_stack[current_ps.depth] = -1 + # If every choice of point to split gave something in the + # (same!) primary orbit, then all children of the first + # stack at this point are equivalent. + j = 0 + for i from 0 <= i < n: + if bitset_check(vertices_to_split[current_ps.depth], i): + j += 1 + if j == subgroup_primary_orbit_size and first_kids_are_same == current_ps.depth+1: + first_kids_are_same = current_ps.depth + # Update group size and backtrack. + mpz_mul_si(subgroup_size, subgroup_size, subgroup_primary_orbit_size) + subgroup_primary_orbit_size = 0 + current_ps.depth -= 1 + break + if not new_vertex: + continue + + if current_kids_are_same > current_ps.depth + 1: + current_kids_are_same = current_ps.depth + 1 + if first_and_current_indicator_same > current_ps.depth: + first_and_current_indicator_same = current_ps.depth + if canonical_label and label_and_current_indicator_same > current_ps.depth: + label_and_current_indicator_same = current_ps.depth + compared_current_and_label_indicators = 0 + + # II. Refine down to a discrete partition, or until + # we leave the part of the tree we are interested in + discrete = 0 + while 1: + i = current_ps.depth + current_indicators[i] = split_point_and_refine(current_ps, vertices_determining_current_stack[i], S, refine_and_return_invariant, cells_to_refine_by) + if first_and_current_indicator_same == i and current_indicators[i] == first_indicators[i]: + first_and_current_indicator_same += 1 + if canonical_label: + if compared_current_and_label_indicators == 0: + if label_indicators[i] == -1: + compared_current_and_label_indicators = -1 + else: + compared_current_and_label_indicators = cmp(current_indicators[i], label_indicators[i]) + if label_and_current_indicator_same == i and compared_current_and_label_indicators == 0: + label_and_current_indicator_same += 1 + if compared_current_and_label_indicators > 0: + label_indicators[i] = current_indicators[i] + if first_and_current_indicator_same < current_ps.depth and (not canonical_label or compared_current_and_label_indicators < 0): + break + if PS_is_discrete(current_ps): + discrete = 1 + break + vertices_determining_current_stack[current_ps.depth] = PS_first_smallest(current_ps, vertices_to_split[current_ps.depth]) + bitset_unset(vertices_have_been_reduced, current_ps.depth) + if not all_children_are_equivalent(current_ps, S): + current_kids_are_same = current_ps.depth + 1 + + # III. Check for automorphisms and labels + automorphism = 0 + backtrack = 1 - discrete + if discrete: + if current_ps.depth == first_and_current_indicator_same: + PS_get_perm_from(current_ps, first_ps, permutation) + if compare_structures(permutation, id_perm, S) == 0: + automorphism = 1 + if not automorphism: + if (not canonical_label) or (compared_current_and_label_indicators < 0): + backtrack = 1 + else: + # if we get here, discrete must be true + update_label = 0 + if (compared_current_and_label_indicators > 0) or (current_ps.depth < label_ps.depth): + update_label = 1 + else: + i = compare_structures(current_ps.entries, label_ps.entries, S) + if i > 0: + update_label = 1 + elif i < 0: + backtrack = 1 + else: + PS_get_perm_from(current_ps, label_ps, permutation) + automorphism = 1 + if update_label: + PS_copy_from_to(current_ps, label_ps) + compared_current_and_label_indicators = 0 + label_meets_current = current_ps.depth + label_and_current_indicator_same = current_ps.depth + label_indicators[current_ps.depth] = -1 + backtrack = 1 + if automorphism: + if index_in_fp_and_mcr < len_of_fp_and_mcr - 1: + index_in_fp_and_mcr += 1 + 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) + if OP_merge_list_perm(orbits_of_subgroup, permutation): # if permutation made orbits coarser + if n*output.num_gens == size_of_generator_array: + # must double its size + size_of_generator_array *= 2 + output.generators = sage_realloc( output.generators, size_of_generator_array ) + if output.generators is NULL: + mem_err = True + break # main loop + j = n*output.num_gens + for i from 0 <= i < n: + output.generators[j + i] = permutation[i] + output.num_gens += 1 + if orbits_of_subgroup.mcr[OP_find(orbits_of_subgroup, minimal_in_primary_orbit)] != minimal_in_primary_orbit: + current_ps.depth = first_meets_current + continue # main loop + if canonical_label: + current_ps.depth = label_meets_current + else: + current_ps.depth = first_meets_current + if bitset_check(vertices_have_been_reduced, current_ps.depth): + bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + continue # main loop + if backtrack: + i = current_ps.depth + current_ps.depth = min(max(first_kids_are_same - 1, label_and_current_indicator_same), current_kids_are_same - 1) + if i == current_kids_are_same: + continue # main loop + if index_in_fp_and_mcr < len_of_fp_and_mcr - 1: + index_in_fp_and_mcr += 1 + bitset_zero(fixed_points_of_generators[index_in_fp_and_mcr]) + bitset_zero(minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + j = current_ps.depth + current_ps.depth = i # just for mcr and fixed functions... + for i from 0 <= i < n: + if PS_is_mcr(current_ps, i): + bitset_set(minimal_cell_reps_of_generators[index_in_fp_and_mcr], i) + if PS_is_fixed(current_ps, i): + bitset_set(fixed_points_of_generators[index_in_fp_and_mcr], i) + current_ps.depth = j + if bitset_check(vertices_have_been_reduced, current_ps.depth): + bitset_and(vertices_to_split[current_ps.depth], vertices_to_split[current_ps.depth], minimal_cell_reps_of_generators[index_in_fp_and_mcr]) + + # End of main loop. + + mpz_init_set(output.order, subgroup_size) + if canonical_label: + for i from 0 <= i < n: + output.relabeling[label_ps.entries[i]] = i + + # Deallocate: + for i from 0 <= i < len_of_fp_and_mcr: + bitset_clear(fixed_points_of_generators[i]) + bitset_clear(minimal_cell_reps_of_generators[i]) + for i from 0 <= i < n: + bitset_clear(vertices_to_split[i]) + bitset_clear(vertices_have_been_reduced) + sage_free(current_indicators) + sage_free(first_indicators) + if canonical_label: + PS_dealloc(label_ps) + sage_free(label_indicators) + sage_free(fixed_points_of_generators) + sage_free(minimal_cell_reps_of_generators) + sage_free(vertices_to_split) + sage_free(permutation) + sage_free(id_perm) + sage_free(cells_to_refine_by) + sage_free(vertices_determining_current_stack) + PS_dealloc(current_ps) + PS_dealloc(first_ps) + OP_dealloc(orbits_of_subgroup) + mpz_clear(subgroup_size) + if not mem_err: + return output + else: + mpz_clear(output.order) + sage_free(output.generators) + if base: + sage_free(output.base) + if canonical_label: + sage_free(output.relabeling) + sage_free(output) + output = NULL + raise MemoryError + + + + + + + + + + + + + + diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/refinement_graphs.pxd --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/refinement_graphs.pxd Mon Aug 11 18:26:07 2008 -0700 @@ -0,0 +1,31 @@ + +#***************************************************************************** +# 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' + +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 + +cdef class GraphStruct: + cdef CGraph G + cdef bint directed + cdef bint use_indicator + cdef int *scratch # length 3n+1 + +cdef int refine_by_degree(PartitionStack *, object, int *, int) +cdef int compare_graphs(int *, int *, object) +cdef bint all_children_are_equivalent(PartitionStack *, object) +cdef inline int degree(PartitionStack *, CGraph, int, int, bint) +cdef inline int sort_by_function(PartitionStack *, int, int *) + + + diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/refinement_graphs.pyx --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/groups/perm_gps/partn_ref/refinement_graphs.pyx Mon Aug 11 18:26:07 2008 -0700 @@ -0,0 +1,699 @@ +""" +Module for graph-theoretic partition backtrack functions. + +DOCTEST: + sage: import sage.groups.perm_gps.partn_ref.refinement_graphs + +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/ +#***************************************************************************** + +def search_tree(G_in, partition, lab=True, dig=False, dict_rep=False, certify=False, + verbosity=0, use_indicator_function=True, sparse=False, + base=False, order=False): + """ + Compute canonical labels and automorphism groups of graphs. + + INPUT: + G_in -- a Sage graph + partition -- a list of lists representing a partition of the vertices + lab -- if True, compute and return the canonical label in addition to the + automorphism group. + dig -- set to True for digraphs and graphs with loops. If True, does not + use optimizations based on Lemma 2.25 in [1] that are valid only for + simple graphs. + dict_rep -- if True, return a dictionary with keys the vertices of the + input graph G_in and values elements of the set the permutation group + acts on. (The point is that graphs are arbitrarily labelled, often + 0..n-1, and permutation groups always act on 1..n. This dictionary + maps vertex labels (such as 0..n-1) to the domain of the permutations.) + certify -- if True, return the permutation from G to its canonical label. + verbosity -- currently ignored + use_indicator_function -- option to turn off indicator function + (True is generally faster) + sparse -- whether to use sparse or dense representation of the graph + (ignored if G is already a CGraph - see sage.graphs.base) + base -- whether to return the first sequence of split vertices (used in + computing the order of the group) + order -- whether to return the order of the automorphism group + + OUTPUT: + Depends on the options. If more than one thing is returned, they are in a + tuple in the following order: + list of generators in list-permutation format -- always + dict -- if dict_rep + graph -- if lab + dict -- if certify + list -- if base + integer -- if order + + DOCTEST: + sage: st = sage.groups.perm_gps.partn_ref.refinement_graphs.search_tree + sage: from sage.graphs.base.dense_graph import DenseGraph + sage: from sage.graphs.base.sparse_graph import SparseGraph + + Graphs on zero vertices: + sage: G = Graph() + sage: st(G, [[]], order=True) + ([], Graph on 0 vertices, 1) + + Graphs on one vertex: + sage: G = Graph(1) + sage: st(G, [[0]], order=True) + ([], Graph on 1 vertex, 1) + + Graphs on two vertices: + sage: G = Graph(2) + sage: st(G, [[0,1]], order=True) + ([[1, 0]], Graph on 2 vertices, 2) + sage: st(G, [[0],[1]], order=True) + ([], Graph on 2 vertices, 1) + sage: G.add_edge(0,1) + sage: st(G, [[0,1]], order=True) + ([[1, 0]], Graph on 2 vertices, 2) + sage: st(G, [[0],[1]], order=True) + ([], Graph on 2 vertices, 1) + + Graphs on three vertices: + sage: G = Graph(3) + sage: st(G, [[0,1,2]], order=True) + ([[0, 2, 1], [1, 0, 2]], Graph on 3 vertices, 6) + sage: st(G, [[0],[1,2]], order=True) + ([[0, 2, 1]], Graph on 3 vertices, 2) + sage: st(G, [[0],[1],[2]], order=True) + ([], Graph on 3 vertices, 1) + sage: G.add_edge(0,1) + sage: st(G, [range(3)], order=True) + ([[1, 0, 2]], Graph on 3 vertices, 2) + sage: st(G, [[0],[1,2]], order=True) + ([], Graph on 3 vertices, 1) + sage: st(G, [[0,1],[2]], order=True) + ([[1, 0, 2]], Graph on 3 vertices, 2) + + The Dodecahedron has automorphism group of size 120: + sage: G = graphs.DodecahedralGraph() + sage: Pi = [range(20)] + sage: st(G, Pi, order=True)[2] + 120 + + The three-cube has automorphism group of size 48: + sage: G = graphs.CubeGraph(3) + sage: G.relabel() + sage: Pi = [G.vertices()] + sage: st(G, Pi, order=True)[2] + 48 + + We obtain the same output using different types of Sage graphs: + sage: G = graphs.DodecahedralGraph() + sage: GD = DenseGraph(20) + sage: GS = SparseGraph(20) + sage: for i,j,_ in G.edge_iterator(): + ... GD.add_arc(i,j); GD.add_arc(j,i) + ... GS.add_arc(i,j); GS.add_arc(j,i) + sage: Pi=[range(20)] + sage: a,b = st(G, Pi) + sage: asp,bsp = st(GS, Pi) + sage: ade,bde = st(GD, Pi) + sage: bsg = Graph(implementation='networkx') + sage: bdg = Graph(implementation='networkx') + sage: for i in range(20): + ... for j in range(20): + ... if bsp.has_arc(i,j): + ... bsg.add_edge(i,j) + ... if bde.has_arc(i,j): + ... bdg.add_edge(i,j) + sage: print a, b.graph6_string() + [[0, 19, 3, 2, 6, 5, 4, 17, 18, 11, 10, 9, 13, 12, 16, 15, 14, 7, 8, 1], [0, 1, 8, 9, 13, 14, 7, 6, 2, 3, 19, 18, 17, 4, 5, 15, 16, 12, 11, 10], [1, 8, 9, 10, 11, 12, 13, 14, 7, 6, 2, 3, 4, 5, 15, 16, 17, 18, 19, 0]] S?[PG__OQ@?_?_?P?CO?_?AE?EC?Ac?@O + sage: a == asp + True + sage: a == ade + True + sage: b == bsg + True + sage: b == bdg + True + + Cubes! + sage: C = graphs.CubeGraph(1) + sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order + 2 + sage: C = graphs.CubeGraph(2) + sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order + 8 + sage: C = graphs.CubeGraph(3) + sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order + 48 + sage: C = graphs.CubeGraph(4) + sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order + 384 + sage: C = graphs.CubeGraph(5) + sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order + 3840 + sage: C = graphs.CubeGraph(6) + sage: gens, order = st(C, [C.vertices()], lab=False, order=True); order + 46080 + + One can also turn off the indicator function (note- this will take longer) + sage: D1 = DiGraph({0:[2],2:[0],1:[1]}, loops=True) + sage: D2 = DiGraph({1:[2],2:[1],0:[0]}, loops=True) + sage: a,b = st(D1, [D1.vertices()], dig=True, use_indicator_function=False) + sage: c,d = st(D2, [D2.vertices()], dig=True, use_indicator_function=False) + sage: b==d + True + + This example is due to Chris Godsil: + sage: HS = graphs.HoffmanSingletonGraph() + sage: clqs = (HS.complement()).cliques() + sage: alqs = [Set(c) for c in clqs if len(c) == 15] + sage: Y = Graph([alqs, lambda s,t: len(s.intersection(t))==0], implementation='networkx') + sage: Y0,Y1 = Y.connected_components_subgraphs() + sage: st(Y0, [Y0.vertices()])[1] == st(Y1, [Y1.vertices()])[1] + True + sage: st(Y0, [Y0.vertices()])[1] == st(HS, [HS.vertices()])[1] + True + sage: st(HS, [HS.vertices()])[1] == st(Y1, [Y1.vertices()])[1] + True + + Certain border cases need to be tested as well: + sage: G = Graph('Fll^G') + sage: a,b,c = st(G, [range(G.num_verts())], order=True); b + Graph on 7 vertices + sage: c + 48 + sage: G = Graph(21) + sage: st(G, [range(G.num_verts())], order=True)[2] == factorial(21) + True + + sage: G = Graph('^????????????????????{??N??@w??FaGa?PCO@CP?AGa?_QO?Q@G?CcA??cc????Bo????{????F_') + sage: perm = {3:15, 15:3} + sage: H = G.relabel(perm, inplace=False) + sage: st(G, [range(G.num_verts())])[1] == st(H, [range(H.num_verts())])[1] + True + + """ + cdef CGraph G + cdef int i, j, n + cdef Integer I + cdef aut_gp_and_can_lab_return *output + cdef int **part + from sage.graphs.graph import GenericGraph, Graph, DiGraph + if isinstance(G_in, GenericGraph): + n = G_in.num_verts() + G_in = G_in.copy() + if G_in.vertices() != range(n): + to = G_in.relabel(return_map=True) + frm = {} + for v in to.iterkeys(): + frm[to[v]] = v + partition = [[to[v] for v in cell] for cell in partition] + else: + to = range(n) + frm = to + if sparse: + G = SparseGraph(n) + else: + G = DenseGraph(n) + if G_in.is_directed(): + for i from 0 <= i < n: + for _,j,_ in G_in.outgoing_edge_iterator(i): + G.add_arc(i,j) + else: + for i from 0 <= i < n: + for _,j,_ in G_in.edge_iterator(i): + if j <= i: + G.add_arc(i,j) + G.add_arc(j,i) + elif isinstance(G_in, CGraph): + G = G_in + n = G.num_verts + to = range(n) + frm = to + else: + raise TypeError("G must be a Sage graph.") + + part = sage_malloc((len(partition)+1) * sizeof(int *)) + if part is NULL: + 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 + + cdef GraphStruct GS = GraphStruct() + GS.G = G + GS.directed = 1 if dig else 0 + GS.use_indicator = 1 if use_indicator_function else 0 + GS.scratch = sage_malloc( (3*G.num_verts + 1) * sizeof(int) ) + if GS.scratch is NULL: + for j from 0 <= j < len(partition): + sage_free(part[j]) + sage_free(part) + raise MemoryError + + output = get_aut_gp_and_can_lab(GS, part, G.num_verts, &all_children_are_equivalent, &refine_by_degree, &compare_graphs, lab, base, order) + sage_free( GS.scratch ) + # prepare output + list_of_gens = [] + for i from 0 <= i < output.num_gens: + list_of_gens.append([output.generators[j+i*G.num_verts] for j from 0 <= j < G.num_verts]) + return_tuple = [list_of_gens] + if dict_rep: + ddd = {} + for v in frm.iterkeys(): + if frm[v] != 0: + ddd[v] = frm[v] + else: + ddd[v] = n + return_tuple.append(ddd) + if lab: + if isinstance(G_in, GenericGraph): + G_C = G_in.copy() + G_C.relabel([output.relabeling[i] for i from 0 <= i < n]) + else: + if isinstance(G, SparseGraph): + G_C = SparseGraph(n) + else: + G_C = DenseGraph(n) + for i from 0 <= i < n: + for j in G.out_neighbors(i): + G_C.add_arc(output.relabeling[i],output.relabeling[j]) + return_tuple.append(G_C) + if certify: + dd = {} + for i from 0 <= i < G.num_verts: + dd[frm[i]] = output.relabeling[i] + return_tuple.append(dd) + if base: + return_tuple.append([output.base[i] for i from 0 <= i < output.base_size]) + if order: + I = Integer() + mpz_set(I.value, output.order) + return_tuple.append(I) + for i from 0 <= i < len(partition): + sage_free(part[i]) + sage_free(part) + mpz_clear(output.order) + sage_free(output.generators) + if base: + sage_free(output.base) + if lab: + sage_free(output.relabeling) + sage_free(output) + if len(return_tuple) == 1: + return return_tuple[0] + else: + return tuple(return_tuple) + +cdef int refine_by_degree(PartitionStack *PS, object S, int *cells_to_refine_by, int ctrb_len): + r""" + Refines the input partition by checking degrees of vertices to the given + cells. + + INPUT: + PS -- a partition stack, whose finest partition is the partition to be + refined. + S -- a graph struct object, which contains scratch space, the graph in + question, and some flags. + cells_to_refine_by -- a list of pointers to cells to check degrees against + in refining the other cells (updated in place) + ctrb_len -- how many cells in cells_to_refine_by + + OUTPUT: + + An integer invariant under the orbits of $S_n$. That is, if $\gamma$ is a + permutation of the vertices, then + $$ I(G, PS, cells_to_refine_by) = I( \gamma(G), \gamma(PS), \gamma(cells_to_refine_by) ) .$$ + + """ + cdef GraphStruct GS = S + cdef CGraph G = GS.G + cdef int current_cell_against = 0 + cdef int current_cell, i, r + cdef int first_largest_subcell + cdef int invariant = 1 + cdef int *degrees = GS.scratch # length 3n+1 + cdef bint necessary_to_split_cell + cdef int against_index + while not PS_is_discrete(PS) and current_cell_against < ctrb_len: + invariant += 1 + current_cell = 0 + while current_cell < PS.degree: + invariant += 50 + i = current_cell + necessary_to_split_cell = 0 + while 1: + degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 0) + if degrees[i-current_cell] != degrees[0]: + necessary_to_split_cell = 1 + i += 1 + if PS.levels[i-1] <= PS.depth: + break + # now, i points to the next cell (before refinement) + if necessary_to_split_cell: + invariant += 10 + first_largest_subcell = sort_by_function(PS, current_cell, degrees) + invariant += first_largest_subcell + against_index = current_cell_against + while against_index < ctrb_len: + if cells_to_refine_by[against_index] == current_cell: + cells_to_refine_by[against_index] = first_largest_subcell + break + against_index += 1 + r = current_cell + while 1: + if r == current_cell or PS.levels[r-1] == PS.depth: + if r != first_largest_subcell: + cells_to_refine_by[ctrb_len] = r + ctrb_len += 1 + r += 1 + if r >= i: + break + invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 0) + invariant += (i - current_cell) + current_cell = i + if GS.directed: + # if we are looking at a digraph, also compute + # the reverse degrees and sort by them + current_cell = 0 + while current_cell < PS.degree: # current_cell is still a valid cell + invariant += 20 + i = current_cell + necessary_to_split_cell = 0 + while 1: + degrees[i-current_cell] = degree(PS, G, i, cells_to_refine_by[current_cell_against], 1) + if degrees[i-current_cell] != degrees[0]: + necessary_to_split_cell = 1 + i += 1 + if PS.levels[i-1] <= PS.depth: + break + # now, i points to the next cell (before refinement) + if necessary_to_split_cell: + invariant += 7 + first_largest_subcell = sort_by_function(PS, current_cell, degrees) + invariant += first_largest_subcell + against_index = current_cell_against + while against_index < ctrb_len: + if cells_to_refine_by[against_index] == current_cell: + cells_to_refine_by[against_index] = first_largest_subcell + break + against_index += 1 + against_index = ctrb_len + r = current_cell + while 1: + if r == current_cell or PS.levels[r-1] == PS.depth: + if r != first_largest_subcell: + cells_to_refine_by[against_index] = r + against_index += 1 + ctrb_len += 1 + r += 1 + if r >= i: + break + invariant += degree(PS, G, i-1, cells_to_refine_by[current_cell_against], 1) + invariant += (i - current_cell) + current_cell = i + current_cell_against += 1 + if GS.use_indicator: + return invariant + else: + return 0 + +cdef int compare_graphs(int *gamma_1, int *gamma_2, object S): + r""" + Compare gamma_1(G) and gamma_2(G). + + 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 + \code{cmp}) function. + + INPUT: + gamma_1, gamma_2 -- list permutations (inverse) + S -- a graph struct object + + """ + 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]): + return 1 + elif G.has_arc_unsafe(gamma_2[i], gamma_2[j]): + return -1 + return 0 + +cdef bint all_children_are_equivalent(PartitionStack *PS, object S): + """ + Return True if every refinement of the current partition results in the + same structure. + + WARNING: + Converse does not hold in general! See Lemma 2.25 of [1] for details. + + INPUT: + PS -- the partition stack to be checked + S -- a graph struct object + """ + cdef GraphStruct GS = S + cdef CGraph G = GS.G + cdef int i, n = PS.degree + cdef bint in_cell = 0 + cdef int nontrivial_cells = 0 + cdef int total_cells = PS_num_cells(PS) + if n <= total_cells + 4: + return 1 + for i from 0 <= i < n-1: + if PS.levels[i] <= PS.depth: + if in_cell: + nontrivial_cells += 1 + in_cell = 0 + else: + in_cell = 1 + if in_cell: + nontrivial_cells += 1 + if n == total_cells + nontrivial_cells: + return 1 + if n == total_cells + nontrivial_cells + 1: + return 1 + return 0 + +cdef inline int degree(PartitionStack *PS, CGraph G, int entry, int cell_index, bint reverse): + """ + Returns the number of edges from the vertex corresponding to entry to + vertices in the cell corresponding to cell_index. + + INPUT: + PS -- the partition stack to be checked + S -- a graph struct object + entry -- the position of the vertex in question in the entries of PS + cell_index -- the starting position of the cell in question in the entries + of PS + reverse -- whether to check for arcs in the other direction + """ + cdef int num_arcs = 0 + entry = PS.entries[entry] + if not reverse: + while 1: + if G.has_arc_unsafe(PS.entries[cell_index], entry): + num_arcs += 1 + if PS.levels[cell_index] > PS.depth: + cell_index += 1 + else: + break + else: + while 1: + if G.has_arc_unsafe(entry, PS.entries[cell_index]): + num_arcs += 1 + if PS.levels[cell_index] > PS.depth: + cell_index += 1 + else: + break + return num_arcs + +cdef inline int sort_by_function(PartitionStack *PS, int start, int *degrees): + """ + A simple counting sort, given the degrees of vertices to a certain cell. + + INPUT: + PS -- the partition stack to be checked + start -- beginning index of the cell to be sorted + degrees -- the values to be sorted by + + """ + cdef int n = PS.degree + cdef int i, j, max, max_location + cdef int *counts = degrees + n, *output = degrees + 2*n + 1 + for i from 0 <= i <= n: + counts[i] = 0 + i = 0 + while PS.levels[i+start] > PS.depth: + counts[degrees[i]] += 1 + i += 1 + counts[degrees[i]] += 1 + # i+start is the right endpoint of the cell now + max = counts[0] + max_location = 0 + for j from 0 < j <= n: + if counts[j] > max: + max = counts[j] + max_location = j + counts[j] += counts[j - 1] + for j from i >= j >= 0: + counts[degrees[j]] -= 1 + output[counts[degrees[j]]] = PS.entries[start+j] + max_location = counts[max_location]+start + for j from 0 <= j <= i: + PS.entries[start+j] = output[j] + j = 1 + while j <= n and counts[j] <= i: + if counts[j] > 0: + PS.levels[start + counts[j] - 1] = PS.depth + PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1) + j += 1 + return max_location + + +def all_labeled_graphs(n): + """ + Returns all labeled graphs on n vertices {0,1,...,n-1}. Used in + classifying isomorphism types (naive approach), and more importantly + in benchmarking the search algorithm. + + EXAMPLE: + sage: from sage.groups.perm_gps.partn_ref.refinement_graphs import all_labeled_graphs + sage: st = sage.groups.perm_gps.partn_ref.refinement_graphs.search_tree + sage: Glist = {} + sage: Giso = {} + sage: for n in [1..5]: + ... Glist[n] = all_labeled_graphs(n) + ... Giso[n] = [] + ... for g in Glist[n]: + ... a, b = st(g, [range(n)]) + ... inn = False + ... for gi in Giso[n]: + ... if b == gi: + ... inn = True + ... if not inn: + ... Giso[n].append(b) + sage: for n in Giso: + ... print n, len(Giso[n]) + 1 1 + 2 2 + 3 4 + 4 11 + 5 34 + + """ + from sage.graphs.graph import Graph + TE = [] + for i in range(n): + for j in range(i): + TE.append((i, j)) + m = len(TE) + Glist= [] + for i in range(2**m): + G = Graph(n) + b = Integer(i).binary() + b = '0'*(m-len(b)) + b + for i in range(m): + if int(b[i]): + G.add_edge(TE[i]) + Glist.append(G) + return Glist + + +def random_tests(t=10.0, n_max=60, perms_per_graph=10): + """ + Tests to make sure that C(gamma(G)) == C(G) for random permutations gamma + and random graphs G. + + INPUT: + t -- run tests for approximately this many seconds + n_max -- test graphs with at most this many vertices + perms_per_graph -- test each graph with this many random permutations + + DISCUSSION: + + Until t seconds have elapsed, this code generates a random graph G on at + most n_max vertices. The density of edges is chosen randomly between 0 + and 1. + + 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. + + TESTS: + sage: import sage.groups.perm_gps.partn_ref.refinement_graphs + sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests() + All passed: ... random tests on ... graphs. + sage: sage.groups.perm_gps.partn_ref.refinement_graphs.random_tests(180.0, 200, 30) # long time + All passed: ... random tests on ... graphs. + + """ + from sage.misc.misc import walltime + from sage.misc.prandom import random, randint + from sage.graphs.graph_generators import GraphGenerators, DiGraphGenerators + from sage.combinat.permutation import Permutations + cdef int i, j, num_tests = 0, num_graphs = 0, passed = 1 + GG = GraphGenerators() + DGG = DiGraphGenerators() + t_0 = walltime() + while walltime(t_0) < t: + p = random() + n = randint(1, n_max) + S = Permutations(n) + + G = GG.RandomGNP(n, p) + H = G.copy() + for i from 0 <= i < perms_per_graph: + G1 = search_tree(G, [G.vertices()])[1] + perm = list(S.random_element()) + perm = [perm[j]-1 for j from 0 <= j < n] + G.relabel(perm) + G2 = search_tree(G, [G.vertices()])[1] + if G1 != G2: + print "FAILURE: graph6-" + print H.graph6_string() + print perm + passed = 0 + + D = DGG.RandomDirectedGNP(n, p) + D.loops(True) + for i from 0 <= i < n: + if random() < p: + D.add_edge(i,i) + E = D.copy() + for i from 0 <= i < perms_per_graph: + D1 = search_tree(D, [D.vertices()], dig=True)[1] + perm = list(S.random_element()) + perm = [perm[j]-1 for j from 0 <= j < n] + D.relabel(perm) + D2 = search_tree(D, [D.vertices()], dig=True)[1] + if D1 != D2: + print "FAILURE: dig6-" + print E.dig6_string() + print perm + passed = 0 + num_tests += 2*perms_per_graph + num_graphs += 2 + if passed: + print "All passed: %d random tests on %d graphs."%(num_tests, num_graphs) + + diff -r 5e2426d277eb sage/misc/bitset.pxi --- a/sage/misc/bitset.pxi Wed Jul 30 06:34:58 2008 -0700 +++ b/sage/misc/bitset.pxi Mon Aug 11 18:26:07 2008 -0700 @@ -16,25 +16,6 @@ # This is a .pxi file so that one can inline functions. Doctests in misc_c. - -include "../ext/stdsage.pxi" - -cdef extern from *: - void *memset(void *, int, size_t) - void *memcpy(void *, void *, size_t) - int memcmp(void *, void *, size_t) - size_t strlen(char *) - - # constant literals - int index_shift "(sizeof(unsigned long)==8 ? 6 : 5)" - unsigned long offset_mask "(sizeof(unsigned long)==8 ? 0x3F : 0x1F)" - -cdef struct bitset_s: - long size - long limbs - unsigned long *bits - -ctypedef bitset_s bitset_t[1] ############################################################################# # Bitset Initalization diff -r 5e2426d277eb sage/misc/bitset_pxd.pxi --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/sage/misc/bitset_pxd.pxi Mon Aug 11 18:26:07 2008 -0700 @@ -0,0 +1,38 @@ +#***************************************************************************** +# Copyright (C) 2008 Robert Bradshaw +# +# Distributed under the terms of the GNU General Public License (GPL) +# +# This code is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# The full text of the GPL is available at: +# +# http://www.gnu.org/licenses/ +#***************************************************************************** + + + +# This is a .pxi file so that one can inline functions. Doctests in misc_c. + +include "../ext/stdsage.pxi" + +cdef extern from *: + void *memset(void *, int, size_t) + void *memcpy(void *, void *, size_t) + int memcmp(void *, void *, size_t) + size_t strlen(char *) + + # constant literals + int index_shift "(sizeof(unsigned long)==8 ? 6 : 5)" + unsigned long offset_mask "(sizeof(unsigned long)==8 ? 0x3F : 0x1F)" + +cdef struct bitset_s: + long size + long limbs + unsigned long *bits + +ctypedef bitset_s bitset_t[1] + diff -r 5e2426d277eb sage/misc/misc_c.pyx --- a/sage/misc/misc_c.pyx Wed Jul 30 06:34:58 2008 -0700 +++ b/sage/misc/misc_c.pyx Mon Aug 11 18:26:07 2008 -0700 @@ -277,6 +277,7 @@ class NonAssociative: # Bitset Testing ############################################################################# +include "bitset_pxd.pxi" include "bitset.pxi" def test_bitset(py_a, py_b, long n): diff -r 5e2426d277eb setup.py --- a/setup.py Wed Jul 30 06:34:58 2008 -0700 +++ b/setup.py Mon Aug 11 18:26:07 2008 -0700 @@ -675,6 +675,12 @@ ext_modules = [ \ Extension('sage.groups.perm_gps.permgroup_element', sources = ['sage/groups/perm_gps/permgroup_element.pyx']), \ + + 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.refinement_graphs', + sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx']), \ Extension('sage.structure.sage_object', sources = ['sage/structure/sage_object.pyx']), \ @@ -1370,6 +1376,7 @@ code = setup(name = 'sage', 'sage.groups.abelian_gps', 'sage.groups.matrix_gps', 'sage.groups.perm_gps', + 'sage.groups.perm_gps.partn_ref', 'sage.interfaces',