Context Navigation

Back to Ticket #3676

Ticket #3676: trac_3676-graph_isom_refactor.patch

File trac_3676-graph_isom_refactor.patch, 88.0 KB (added by rlm, 20 months ago)
Apply only this patch (please do not delete the others!)

sage/graphs/graph_generators.py

# HG changeset patch
# User Robert L. Miller <rlm@rlmiller.org>
# Date 1216425781 25200
# Node ID 074731a2698613667ad1b649a96a8446a892ffd6
# Parent  91af4c3f6b92316e4117af342d82be33cf5d949b
Refactored graph isomorphism code
Split up bitset.pxi to allow use in pxd files.

diff -r 5e2426d277eb sage/graphs/graph_generators.py

a	b
2991	2991	import networkx
2992	2992	return graph.DiGraph(networkx.gnc_graph(n, seed))
2993	2993
	2994	def RandomDirectedGNP(self, n, p):
	2995	r"""
	2996	Returns a random digraph on $n$ nodes. Each edge is inserted
	2997	independently with probability $p$.
	2998
	2999	REFERENCES:
	3000	[1] P. Erdos and A. Renyi, On Random Graphs, Publ. Math. 6, 290 (1959).
	3001	[2] E. N. Gilbert, Random Graphs, Ann. Math. Stat., 30, 1141 (1959).
	3002
	3003	PLOTTING:
	3004	When plotting, this graph will use the default spring-layout
	3005	algorithm, unless a position dictionary is specified.
	3006
	3007	EXAMPLE:
	3008	sage: digraphs.RandomDirectedGNP(10, .2).num_verts()
	3009	10
	3010
	3011	"""
	3012	from random import random
	3013	D = graph.DiGraph(n)
	3014	for i in xrange(n):
	3015	for j in xrange(i):
	3016	if random() < p:
	3017	D.add_edge(i,j)
	3018	for j in xrange(i+1,n):
	3019	if random() < p:
	3020	D.add_edge(i,j)
	3021	return D
	3022
2994	3023	def RandomDirectedGNR(self, n, p, seed=None):
2995	3024	"""
2996	3025	Returns a random GNR (growing network with redirection) digraph with n

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/init.py

diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/__init__.py

a	b
	1

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd

diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pxd

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

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx

diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx

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

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/refinement_graphs.pxd

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

(a) /dev/null vs. (b) b/sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

diff -r 5e2426d277eb sage/groups/perm_gps/partn_ref/refinement_graphs.pyx

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

sage/misc/bitset.pxi

diff -r 5e2426d277eb sage/misc/bitset.pxi

a	b
16	16
17	17
18	18	# This is a .pxi file so that one can inline functions. Doctests in misc_c.
19
20		~~include "../ext/stdsage.pxi"~~
21
22		~~cdef extern from *:~~
23		~~void memset(void , int, size_t)~~
24		~~void memcpy(void , void *, size_t)~~
25		~~int memcmp(void , void , size_t)~~
26		~~size_t strlen(char *)~~
27
28		~~# constant literals~~
29		~~int index_shift "(sizeof(unsigned long)==8 ? 6 : 5)"~~
30		~~unsigned long offset_mask "(sizeof(unsigned long)==8 ? 0x3F : 0x1F)"~~
31
32		~~cdef struct bitset_s:~~
33		~~long size~~
34		~~long limbs~~
35		~~unsigned long *bits~~
36
37		~~ctypedef bitset_s bitset_t[1]~~
38	19
39	20	#############################################################################
40	21	# Bitset Initalization

(a) /dev/null vs. (b) b/sage/misc/bitset_pxd.pxi

diff -r 5e2426d277eb sage/misc/bitset_pxd.pxi

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

sage/misc/misc_c.pyx

diff -r 5e2426d277eb sage/misc/misc_c.pyx

a	b
277	277	# Bitset Testing
278	278	#############################################################################
279	279
	280	include "bitset_pxd.pxi"
280	281	include "bitset.pxi"
281	282
282	283	def test_bitset(py_a, py_b, long n):

setup.py

diff -r 5e2426d277eb setup.py

a	b
675	675
676	676	Extension('sage.groups.perm_gps.permgroup_element',
677	677	sources = ['sage/groups/perm_gps/permgroup_element.pyx']), \
	678
	679	Extension('sage.groups.perm_gps.partn_ref.automorphism_group_canonical_label',
	680	sources = ['sage/groups/perm_gps/partn_ref/automorphism_group_canonical_label.pyx']), \
	681
	682	Extension('sage.groups.perm_gps.partn_ref.refinement_graphs',
	683	sources = ['sage/groups/perm_gps/partn_ref/refinement_graphs.pyx']), \
678	684
679	685	Extension('sage.structure.sage_object',
680	686	sources = ['sage/structure/sage_object.pyx']), \
…	…
1370	1376	'sage.groups.abelian_gps',
1371	1377	'sage.groups.matrix_gps',
1372	1378	'sage.groups.perm_gps',
	1379	'sage.groups.perm_gps.partn_ref',
1373	1380
1374	1381	'sage.interfaces',
1375	1382

Download in other formats:

Original Format