Context Navigation

source: sage/graphs/graph_isom.pyx @ 7119:c5830239ad91

Visit:

Revision 7119:c5830239ad91, 86.7 KB checked in by Robert L Miller <rlm@…>, 6 years ago (diff)
more legibility in graph_isom

Line
1	"""
2	N.I.C.E. - Nice (as in open source) Isomorphism Check Engine
3
4	Automorphism group computation and isomorphism checking for graphs.
5
6	This is an open source implementation of Brendan McKay's algorithm for graph
7	automorphism and isomorphism. McKay released a C version of his algorithm,
8	named nauty (No AUTomorphisms, Yes?) under a license that is not GPL
9	compatible. Although the program is open source, reading the source disallows
10	anyone from recreating anything similar and releasing it under the GPL. Also,
11	many people have complained that the code is difficult to understand. The
12	first main goal of NICE was to produce a genuinely open graph isomorphism
13	program, which has been accomplished. The second goal is for this code to be
14	understandable, so that computed results can be trusted and further derived
15	work will be possible.
16
17	To determine the isomorphism type of a graph, it is convenient to define a
18	canonical label for each isomorphism class- essentially an equivalence class
19	representative. Loosely (albeit incorrectly), the canonical label is defined
20	by enumerating all labeled graphs, then picking the maximal one in each
21	isomorphism class. The NICE algorithm is essentially a backtrack search. It
22	searches through the rooted tree of partition nests (where each partition is
23	equitable) for implicit and explicit automorphisms, and uses this information
24	to eliminate large parts of the tree from further searching. Since the leaves
25	of the search tree are all discrete ordered partitions, each one of these
26	corresponds to an ordering of the vertices of the graph, i.e. another member
27	of the isomorphism class. Once the algorithm has finished searching the tree,
28	it will know which leaf corresponds to the canonical label. In the process,
29	generators for the automorphism group are also produced.
30
31	AUTHORS:
32	Robert L. Miller -- (2007-03-20) initial version
33	Tom Boothby -- (2007-03-20) help with indicator function
34	Robert L. Miller -- (2007-04-07--30) optimizations
35	(2007-07-07--14) PartitionStack and OrbitPartition
36	Tom Boothby -- (2007-07-14) datastructure advice
37	Robert L. Miller -- (2007-07-16--20) bug fixes
38
39	REFERENCE:
40	[1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
41	Vol. 30 (1981), pp. 45-87.
42
43	NOTE:
44	Often we assume that G is a graph on vertices {0,1,...,n-1}.
45	"""
46
47	#*****************************************************************************
48	# Copyright (C) 2006 - 2007 Robert L. Miller <rlmillster@gmail.com>
49	#
50	# Distributed under the terms of the GNU General Public License (GPL)
51	# http://www.gnu.org/licenses/
52	#*****************************************************************************
53
54	include '../ext/cdefs.pxi'
55	include '../ext/python_mem.pxi'
56	include '../ext/stdsage.pxi'
57
58	from sage.graphs.graph import Graph, DiGraph
59	from sage.misc.misc import cputime
60	from sage.rings.integer import Integer
61
62	cdef class OrbitPartition:
63	"""
64	TODO: documentation
65
66	EXAMPLES:
67	sage: from sage.graphs.graph_isom import OrbitPartition
68	sage: K = OrbitPartition(20)
69	sage: K.find(7)
70	7
71	sage: K.union_find(7, 12)
72	sage: K.find(12)
73	7
74	sage: J = OrbitPartition(20)
75	sage: J.is_finer_than(K, 20)
76	True
77	sage: K.is_finer_than(J, 20)
78	False
79
80	sage: from sage.graphs.graph_isom import OrbitPartition
81	sage: Theta1 = OrbitPartition(10)
82	sage: Theta2 = OrbitPartition(10)
83	sage: Theta1.union_find(0,1)
84	sage: Theta1.union_find(2,3)
85	sage: Theta1.union_find(3,4)
86	sage: Theta1.union_find(5,6)
87	sage: Theta1.union_find(8,9)
88	sage: Theta2.vee_with(Theta1, 10)
89	sage: for i in range(10):
90	... print i, Theta2.find(i)
91	0 0
92	1 0
93	2 2
94	3 2
95	4 2
96	5 5
97	6 5
98	7 7
99	8 8
100	9 8
101
102	"""
103
104	def __new__(self, int n):
105	cdef int k
106	self.elements = <int > sage_malloc( n sizeof(int) )
107	if not self.elements:
108	raise MemoryError("Error allocating memory.")
109	self.sizes = <int > sage_malloc( n sizeof(int) )
110	if not self.sizes:
111	sage_free(self.elements)
112	raise MemoryError("Error allocating memory.")
113	for k from 0 <= k < n:
114	self.elements[k] = -1
115	self.sizes[k] = 1
116
117	def __dealloc__(self):
118	sage_free(self.elements)
119	sage_free(self.sizes)
120
121	def find(self, x):
122	return self._find(x)
123
124	cdef int _find(self, int x):
125	if self.elements[x] == -1:
126	return x
127	self.elements[x] = self._find(self.elements[x])
128	return self.elements[x]
129
130	def union_find(self, a, b):
131	self._union_find(a, b)
132
133	cdef void _union_find(self, int a, int b):
134	cdef int aRoot, bRoot
135	aRoot = self._find(a)
136	bRoot = self._find(b)
137	self._union_roots(aRoot, bRoot)
138
139	def union_roots(self, a, b):
140	self._union_roots(a, b)
141
142	cdef void _union_roots(self, int a, int b):
143	if a < b:
144	self.elements[b] = a
145	self.sizes[b] += self.sizes[a]
146	elif a > b:
147	self.elements[a] = b
148	self.sizes[a] += self.sizes[b]
149
150	def is_finer_than(self, other, n):
151	return self._is_finer_than(other, n) == 1
152
153	cdef int _is_finer_than(self, OrbitPartition other, int n):
154	cdef int i
155	for i from 0 <= i < n:
156	if self.elements[i] != -1 and other.find(self.find(i)) != other.find(i):
157	return 0
158	return 1
159
160	def vee_with(self, other, n):
161	self._vee_with(other, n)
162
163	cdef void _vee_with(self, OrbitPartition other, int n):
164	cdef int i
165	for i from 0 <= i < n:
166	if self.elements[i] == -1:
167	self._union_roots(i, self.find(other.find(i)))
168
169	cdef int _is_min_cell_rep(self, int i):
170	if self.elements[i] == -1:
171	return 1
172	return 0
173
174	cdef int _is_fixed(self, int i):
175	if self.elements[i] == -1 and self.sizes[i] == 1:
176	return 1
177	return 0
178
179	cdef OrbitPartition _orbit_partition_from_list_perm(int *gamma, int n):
180	cdef int i
181	cdef OrbitPartition O
182	O = OrbitPartition(n)
183	for i from 0 <= i < n:
184	if i != gamma[i]:
185	O._union_find(i, gamma[i])
186	return O
187
188	cdef class PartitionStack:
189	"""
190	TODO: documentation
191
192	EXAMPLES:
193
194	sage: from sage.graphs.graph_isom import PartitionStack
195	sage: P = PartitionStack([range(9, -1, -1)])
196	sage: P.sort_by_function(0, [2,1,2,1,2,1,3,4,2,1], 1, 10)
197	0
198	sage: P.sort_by_function(0, [2,1,2,1], 2, 10)
199	0
200	sage: P.sort_by_function(4, [2,1,2,1], 3, 10)
201	4
202	sage: P.sort_by_function(0, [0,1], 4, 10)
203	0
204	sage: P.sort_by_function(2, [1,0], 5, 10)
205	2
206	sage: P.sort_by_function(4, [1,0], 6, 10)
207	4
208	sage: P.sort_by_function(6, [1,0], 7, 10)
209	6
210	sage: P
211	({5,9,7,1,6,2,8,0,4,3})
212	({5,9,7,1},{6,2,8,0},{4},{3})
213	({5,9},{7,1},{6,2,8,0},{4},{3})
214	({5,9},{7,1},{6,2},{8,0},{4},{3})
215	({5},{9},{7,1},{6,2},{8,0},{4},{3})
216	({5},{9},{7},{1},{6,2},{8,0},{4},{3})
217	({5},{9},{7},{1},{6},{2},{8,0},{4},{3})
218	({5},{9},{7},{1},{6},{2},{8},{0},{4},{3})
219	({5},{9},{7},{1},{6},{2},{8},{0},{4},{3})
220	({5},{9},{7},{1},{6},{2},{8},{0},{4},{3})
221	sage: P.is_discrete(7)
222	1
223	sage: P.is_discrete(6)
224	0
225
226	sage: M = graphs.PetersenGraph().am()
227	sage: MM = []
228	sage: for i in range(10):
229	... MM.append([])
230	... for j in range(10):
231	... MM[i].append(M[i][j])
232	sage: P = PartitionStack(10)
233	sage: P.split_vertex(0, 1)
234	sage: P.refine_by_square_matrix(MM, 1, [0], 10, 0)
235	sage: P
236	({0,2,3,6,7,8,9,1,4,5})
237	({0},{2,3,6,7,8,9},{1,4,5})
238	({0},{2,3,6,7,8,9},{1,4,5})
239	({0},{2,3,6,7,8,9},{1,4,5})
240	({0},{2,3,6,7,8,9},{1,4,5})
241	({0},{2,3,6,7,8,9},{1,4,5})
242	({0},{2,3,6,7,8,9},{1,4,5})
243	({0},{2,3,6,7,8,9},{1,4,5})
244	({0},{2,3,6,7,8,9},{1,4,5})
245	({0},{2,3,6,7,8,9},{1,4,5})
246	sage: P.split_vertex(1, 2)
247	sage: P.refine_by_square_matrix(MM, 2, [7], 10, 0)
248	sage: P
249	({0,3,7,8,9,2,6,1,4,5})
250	({0},{3,7,8,9,2,6},{1,4,5})
251	({0},{3,7,8,9},{2,6},{1},{4,5})
252	({0},{3,7,8,9},{2,6},{1},{4,5})
253	({0},{3,7,8,9},{2,6},{1},{4,5})
254	({0},{3,7,8,9},{2,6},{1},{4,5})
255	({0},{3,7,8,9},{2,6},{1},{4,5})
256	({0},{3,7,8,9},{2,6},{1},{4,5})
257	({0},{3,7,8,9},{2,6},{1},{4,5})
258	({0},{3,7,8,9},{2,6},{1},{4,5})
259
260
261	"""
262	def __new__(self, data):
263	cdef int j, k, n
264	cdef PartitionStack _data
265	try:
266	n = int(data)
267	self.entries = <int > sage_malloc( n sizeof(int) )
268	if not self.entries:
269	raise MemoryError("Error allocating memory.")
270	self.levels = <int > sage_malloc( n sizeof(int) )
271	if not self.levels:
272	sage_free(self.entries)
273	raise MemoryError("Error allocating memory.")
274	for k from 0 <= k < n-1:
275	self.entries[k] = k
276	self.levels[k] = n
277	self.entries[n-1] = n-1
278	self.levels[n-1] = -1
279	except:
280	if isinstance(data, list):
281	n = sum([len(datum) for datum in data])
282	self.entries = <int > sage_malloc( n sizeof(int) )
283	if not self.entries:
284	raise MemoryError("Error allocating memory.")
285	self.levels = <int > sage_malloc( n sizeof(int) )
286	if not self.levels:
287	sage_free(self.entries)
288	raise MemoryError("Error allocating memory.")
289	j = 0
290	k = 0
291	for cell in data:
292	for entry in cell:
293	self.entries[j] = entry
294	self.levels[j] = n
295	j += 1
296	self.levels[j-1] = 0
297	self._percolate(k, j-1)
298	k = j
299	self.levels[j-1] = -1
300	elif isinstance(data, PartitionStack):
301	_data = data
302	j = 0
303	while _data.levels[j] != -1: j += 1
304	n = j + 1
305	self.entries = <int > sage_malloc( n sizeof(int) )
306	if not self.entries:
307	raise MemoryError("Error allocating memory.")
308	self.levels = <int > sage_malloc( n sizeof(int) )
309	if not self.levels:
310	sage_free(self.entries)
311	raise MemoryError("Error allocating memory.")
312	for k from 0 <= k < n:
313	self.entries[k] = _data.entries[k]
314	self.levels[k] = _data.levels[k]
315	else:
316	raise ValueError("Input must be an int, a list of lists, or a PartitionStack.")
317
318	def __dealloc__(self):
319	sage_free(self.entries)
320	sage_free(self.levels)
321
322	def __repr__(self):
323	k = 0
324	s = ''
325	while k == 0 or self.levels[k-1] != -1:
326	s += '({'
327	i = 0
328	while i == 0 or self.levels[i-1] != -1:
329	s += str(self.entries[i])
330	if self.levels[i] <= k:
331	s += '},{'
332	else:
333	s += ','
334	i += 1
335	s = s[:-2] + ')\n'
336	k += 1
337	return s
338
339	def is_discrete(self, k):
340	return self._is_discrete(k)
341
342	cdef int _is_discrete(self, int k):
343	cdef int i = 0
344	while True:
345	if self.levels[i] > k:
346	return 0
347	if self.levels[i] == -1: break
348	i += 1
349	return 1
350
351	def num_cells(self, k):
352	return self._num_cells(k)
353
354	cdef int _num_cells(self, int k):
355	cdef int i = 0, j = 1
356	while self.levels[i] != -1:
357	#for i from 0 <= i < n-1:
358	if self.levels[i] <= k:
359	j += 1
360	i += 1
361	return j
362
363	def sat_225(self, k, n):
364	return self._sat_225(k, n) == 1
365
366	cdef int _sat_225(self, int k, int n):
367	cdef int i, in_cell = 0
368	cdef int nontrivial_cells = 0
369	cdef int total_cells = self._num_cells(k)
370	if n <= total_cells + 4:
371	return 1
372	for i from 0 <= i < n-1:
373	if self.levels[i] <= k:
374	if in_cell:
375	nontrivial_cells += 1
376	in_cell = 0
377	else:
378	in_cell = 1
379	if in_cell:
380	nontrivial_cells += 1
381	if n == total_cells + nontrivial_cells:
382	return 1
383	if n == total_cells + nontrivial_cells + 1:
384	return 1
385	return 0
386
387	cdef int _is_min_cell_rep(self, int i, int k):
388	return i == 0 or self.levels[i-1] <= k
389
390	cdef int _is_fixed(self, int i, int k):
391	"""
392	Assuming you already know it is a minimum cell representative.
393	"""
394	return self.levels[i] <= k
395
396	def split_vertex(self, v, k):
397	"""
398	Splits the cell in self(k) containing v, putting new cells in place
399	in self(k).
400	"""
401	self._split_vertex(v, k)
402
403	cdef int _split_vertex(self, int v, int k):
404	cdef int i = 0, j
405	while self.entries[i] != v:
406	i += 1
407	j = i
408	while self.levels[i] > k:
409	i += 1
410	if j == 0 or self.levels[j-1] <= k:
411	self._percolate(j+1, i)
412	else:
413	while j != 0 and self.levels[j-1] > k:
414	self.entries[j] = self.entries[j-1]
415	j -= 1
416	self.entries[j] = v
417	self.levels[j] = k
418	return j
419
420	def percolate(self, start, end):
421	self._percolate(start, end)
422
423	cdef void _percolate(self, int start, int end):
424	cdef int i, temp
425	for i from end >= i > start:
426	if self.entries[i] < self.entries[i-1]:
427	temp = self.entries[i]
428	self.entries[i] = self.entries[i-1]
429	self.entries[i-1] = temp
430
431	def sort_by_function(self, start, degrees, k, n):
432	cdef int i
433	cdef int degs = <int > sage_malloc( ( 3 * n + 1 ) * sizeof(int) )
434	if not degs:
435	raise MemoryError("Couldn't allocate...")
436	for i from 0 <= i < len(degrees):
437	degs[i] = degrees[i]
438	return self._sort_by_function(start, degs, k, n)
439	sage_free(degs)
440
441	cdef int _sort_by_function(self, int start, int *degrees, int k, int n):
442	cdef int i, j, m = 2*n, max, max_location
443	cdef int counts = degrees + n, output = degrees + 2*n + 1
444	# print '\|'.join(['%02d'%self.entries[iii] for iii in range(n)])
445	# print '\|'.join(['%02d'%self.levels[iii] for iii in range(n)])
446	# print '\|'.join(['%02d'%degrees[iii] for iii in range(n)])
447	# print '\|'.join(['%02d'%counts[iii] for iii in range(n)])
448	# print '\|'.join(['%02d'%output[iii] for iii in range(n)])
449
450	for i from 0 <= i <= n:
451	counts[i] = 0
452	i = 0
453	while self.levels[i+start] > k:
454	counts[degrees[i]] += 1
455	i += 1
456	counts[degrees[i]] += 1
457
458	# i+start is the right endpoint of the cell now
459	max = counts[0]
460	max_location = 0
461	for j from 0 < j <= n:
462	if counts[j] > max:
463	max = counts[j]
464	max_location = j
465	counts[j] += counts[j - 1]
466
467	for j from i >= j >= 0:
468	counts[degrees[j]] -= 1
469	output[counts[degrees[j]]] = self.entries[start+j]
470
471	max_location = counts[max_location]+start
472
473	for j from 0 <= j <= i:
474	self.entries[start+j] = output[j]
475
476	j = 1
477	while j <= n and counts[j] <= i:
478	if counts[j] > 0:
479	self.levels[start + counts[j] - 1] = k
480	self._percolate(start + counts[j-1], start + counts[j] - 1)
481	j += 1
482
483	return max_location
484
485	def clear(self, k):
486	self._clear(k)
487
488	cdef void _clear(self, int k):
489	cdef int i = 0, j = 0
490	while self.levels[i] != -1:
491	if self.levels[i] >= k:
492	self.levels[i] += 1
493	if self.levels[i] < k:
494	self._percolate(j, i)
495	j = i + 1
496	i+=1
497
498	def refine_by_square_matrix(self, G_matrix, k, alpha, n, dig):
499	cdef int *_alpha, i, j
500	cdef int **G
501	_alpha = <int > sage_malloc( ( 4 n + 1 )* sizeof(int) )
502	if not _alpha:
503	raise MemoryError("Memory!")
504	G = <int *> sage_malloc( n sizeof(int*) )
505	if not G:
506	sage_free(_alpha)
507	raise MemoryError("Memory!")
508	for i from 0 <= i < n:
509	G[i] = <int > sage_malloc( n sizeof(int) )
510	if not G[i]:
511	for j from 0 <= j < i:
512	sage_free(G[j])
513	sage_free(G)
514	sage_free(_alpha)
515	raise MemoryError("Memory!")
516	for i from 0 <= i < n:
517	for j from 0 <= j < n:
518	G[i][j] = G_matrix[i][j]
519	# note -- one could memcopy or memmove for the above loop,
520	# but this function is only a python interface to a C function
521	# that gets called by other functions which have already allocated
522	# memory
523	for i from 0 <= i < len(alpha):
524	_alpha[i] = alpha[i]
525	_alpha[len(alpha)] = -1
526	self._refine_by_square_matrix(k, _alpha, n, G, dig)
527	sage_free(_alpha)
528	for i from 0 <= i < n:
529	sage_free(G[i])
530	sage_free(G)
531
532	cdef int _refine_by_square_matrix(self, int k, int alpha, int n, int *G, int dig):
533	cdef int m = 0, j # - m iterates through alpha, the indicator cells
534	# - j iterates through the cells of the partition
535	cdef int i, t, s, r # local variables:
536	# - s plays a double role: outer role indicates whether
537	# splitting the cell is necessary, inner role is as an index
538	# for augmenting _alpha
539	# - i, r iterators
540	# - t: holds the first largest subcell from sort function
541	cdef int invariant = 1
542	# as described in [1], an indicator function Lambda(G, Pi, nu) is
543	# needed to differentiate nonisomorphic branches on the search
544	# tree. The condition is simply that this invariant not depend
545	# on a simultaneous relabeling of the graph G, the root partition
546	# Pi, and the partition nest nu. Since the function will execute
547	# exactly the same way regardless of the labelling, anything that
548	# does not depend on self.entries goes... at least, anything cheap
549	cdef int degrees = alpha + n # alpha assumed to be length 4n + 1 for
550	# extra scratch space
551	while not self._is_discrete(k) and alpha[m] != -1:
552	invariant += 1
553	j = 0
554	while j < n: # j still points at a valid cell
555	invariant += 50
556	# print ' '
557	# print '\|'.join(['%02d'%self.entries[iii] for iii in range(n)])
558	# print '\|'.join(['%02d'%self.levels[iii] for iii in range(n)])
559	# print '\|'.join(['%02d'%alpha[iii] for iii in range(n)])
560	# print '\|'.join(['%02d'%degrees[iii] for iii in range(n)])
561	# print 'j =', j
562	# print 'm =', m
563	i = j; s = 0
564	while True:
565	degrees[i-j] = self._degree_square_matrix(G, i, alpha[m], k)
566	if degrees[i-j] != degrees[0]: s = 1
567	i += 1
568	if self.levels[i-1] <= k: break
569	# print '\|'.join(['%02d'%degrees[iii] for iii in range(n)])
570	# now: j points to this cell,
571	# i points to the next cell (before refinement)
572	if s:
573	invariant += 10
574	t = self._sort_by_function(j, degrees, k, n)
575	# t now points to the first largest subcell
576	invariant += t + degrees[i - j - 1]
577	s = m
578	while alpha[s] != -1:
579	if alpha[s] == j: alpha[s] = t
580	s += 1
581	r = j
582	while True:
583	if r == 0 or self.levels[r-1] == k:
584	if r != t:
585	alpha[s] = r
586	s += 1
587	r += 1
588	if r >= i: break
589	alpha[s] = -1
590	while self.levels[j] > k:
591	j += 1
592	j += 1
593	invariant += (i - j)
594	else: j = i
595	if not dig: m += 1; continue
596	# if we are looking at a digraph, also compute
597	# the reverse degrees and sort by them
598	j = 0
599	while j < n: # j still points at a valid cell
600	invariant += 20
601	# print ' '
602	# print '\|'.join(['%02d'%self.entries[iii] for iii in range(n)])
603	# print '\|'.join(['%02d'%self.levels[iii] for iii in range(n)])
604	# print '\|'.join(['%02d'%alpha[iii] for iii in range(n)])
605	# print '\|'.join(['%02d'%degrees[iii] for iii in range(n)])
606	# print 'j =', j
607	# print 'm =', m
608	i = j; s = 0
609	while True:
610	degrees[i-j] = self._degree_inv_square_matrix(G, i, alpha[m], k)
611	if degrees[i-j] != degrees[0]: s = 1
612	i += 1
613	if self.levels[i-1] <= k: break
614	# now: j points to this cell,
615	# i points to the next cell (before refinement)
616	if s:
617	invariant += 7
618	t = self._sort_by_function(j, degrees, k, n)
619	# t now points to the first largest subcell
620	invariant += t + degrees[i - j - 1]
621	s = m
622	while alpha[s] != -1:
623	if alpha[s] == j: alpha[s] = t
624	s += 1
625	r = j
626	while True:
627	if r == 0 or self.levels[r-1] == k:
628	if r != t:
629	alpha[s] = r
630	s += 1
631	r += 1
632	if r >= i: break
633	alpha[s] = -1
634	while self.levels[j] > k:
635	j += 1
636	j += 1
637	invariant += (i - j)
638	else: j = i
639	m += 1
640	return invariant
641
642	def degree_square_matrix(self, G, v, W, k):
643	cdef int i, j, n = len(G)
644	cdef int GG = <int > sage_malloc( n * sizeof(int*) )
645	if not GG:
646	raise MemoryError("Memory!")
647	for i from 0 <= i < n:
648	GG[i] = <int > sage_malloc( n sizeof(int) )
649	if not GG[i]:
650	for j from 0 <= j < i:
651	sage_free(GG[j])
652	sage_free(GG)
653	raise MemoryError("Memory!")
654	for i from 0 <= i < n:
655	for j from 0 <= j < n:
656	GG[i][j] = G[i][j]
657	j = self._degree_square_matrix(GG, v, W, k)
658	for i from 0 <= i < n:
659	sage_free(GG[i])
660	sage_free(GG)
661	return j
662
663	cdef int _degree_square_matrix(self, int** G, int v, int W, int k):
664	"""
665	G is a square matrix, and W points to the beginning of a cell in the
666	k-th part of the stack.
667	"""
668	cdef int i = 0
669	v = self.entries[v]
670	while True:
671	if G[self.entries[W]][v]:
672	i += 1
673	if self.levels[W] > k: W += 1
674	else: break
675	return i
676
677	cdef int _degree_inv_square_matrix(self, int** G, int v, int W, int k):
678	"""
679	G is a square matrix, and W points to the beginning of a cell in the
680	k-th part of the stack.
681	"""
682	cdef int i = 0
683	v = self.entries[v]
684	while True:
685	if G[v][self.entries[W]]:
686	i += 1
687	if self.levels[W] > k: W += 1
688	else: break
689	return i
690
691	cdef int _first_smallest_nontrivial(self, int *W, int k, int n):
692	cdef int i = 0, j = 0, location = 0, min = n
693	while True:
694	W[i] = 0
695	if self.levels[i] <= k:
696	if i != j and n > i - j + 1:
697	n = i - j + 1
698	location = j
699	j = i + 1
700	if self.levels[i] == -1: break
701	i += 1
702	# location now points to the beginning of the first, smallest,
703	# nontrivial cell
704	while True:
705	if min > self.entries[location]:
706	min = self.entries[location]
707	W[self.entries[location]] = 1
708	if self.levels[location] <= k: break
709	location += 1
710	return min
711
712	cdef void _get_permutation_from(self, PartitionStack zeta, int *gamma):
713	cdef int i = 0
714
715	while True:
716	gamma[zeta.entries[i]] = self.entries[i]
717	i += 1
718	if self.levels[i-1] == -1: break
719
720	# (TODO)
721	# Important note: the enumeration should be kept abstract, and only comparison
722	# functions should be written. This takes up too much memory and time. Simply
723	# iterate starting with the most significant digit in the matrix, and return
724	# as soon as a contradiction is encountered.
725
726	cdef _enumerate_graph_from_discrete(self, int **G, int n):
727	cdef int i, j
728	enumeration = Integer(0)
729	for i from 0 <= i < n:
730	for j from 0 <= j < n:
731	if G[i][j]:
732	enumeration += Integer(2)*((n-(self.entries[i]+1))n + n-(self.entries[j]+1))
733	return enumeration
734
735	cdef _enumerate_graph_with_permutation(int *G, int n, int gamma):
736	cdef int i, j
737	enumeration = Integer(0)
738	for i from 0 <= i < n:
739	for j from 0 <= j < n:
740	if G[i][j]:
741	enumeration += Integer(2)*((n-(gamma[i]+1))n + n-(gamma[j]+1))
742	return enumeration
743
744	cdef _enumerate_graph(int **G, int n):
745	cdef int i, j # enumeration = 0
746	enumeration = Integer(0)
747	for i from 0 <= i < n:
748	for j from 0 <= j < n:
749	if G[i][j]:
750	enumeration += Integer(2)*((n-(i+1))n + n-(j+1))
751	return enumeration
752
753	def _term_pnest_graph(G, PartitionStack nu):
754	"""
755	BDM's G(nu): returns the graph G, relabeled in the order found in
756	nu[m], where m is the first index corresponding to a discrete partition.
757	Assumes nu is a terminal partition nest in T(G, Pi).
758	"""
759	cdef int i, n
760	n = G.order()
761	d = {}
762	for i from 0 <= i < n:
763	d[nu.entries[i]] = i
764	H = G.copy()
765	H.relabel(d)
766	return H
767
768	def search_tree(G, Pi, lab=True, dig=False, dict=False, certify=False, verbosity=0):
769	"""
770	Assumes that the vertex set of G is {0,1,...,n-1} for some n.
771
772	Note that this conflicts with the SymmetricGroup we are using to represent
773	automorphisms. The solution is to let the group act on the set
774	{1,2,...,n}, under the assumption n = 0.
775
776	INPUT:
777	lab-- if True, return the canonical label in addition to the auto-
778	morphism group.
779	dig-- if True, does not use Lemma 2.25 in [1], and the algorithm is
780	valid for digraphs and graphs with loops.
781	dict-- if True, explain which vertices are which elements of the set
782	{1,2,...,n} in the representation of the automorphism group.
783	certify-- if True, return the automorphism between G and its canonical
784	label. Forces lab=True.
785	verbosity-- 0 - print nothing
786	1 - display state trace
787	2 - with timings
788	3 - display partition nests
789	4 - display orbit partition
790
791	STATE DIAGRAM:
792	sage: SD = DiGraph( { 1:[18,2], 2:[5,3], 3:[4,6], 4:[7,2], 5:[4], 6:[13,12], 7:[18,8,10], 8:[6,9,10], 9:[6], 10:[11,13], 11:[12], 12:[13], 13:[17,14], 14:[16,15], 15:[2], 16:[13], 17:[15,13], 18:[13] } )
793	sage: SD.set_edge_label(1, 18, 'discrete')
794	sage: SD.set_edge_label(4, 7, 'discrete')
795	sage: SD.set_edge_label(2, 5, 'h = 0')
796	sage: SD.set_edge_label(7, 18, 'h = 0')
797	sage: SD.set_edge_label(7, 10, 'aut')
798	sage: SD.set_edge_label(8, 10, 'aut')
799	sage: SD.set_edge_label(8, 9, 'label')
800	sage: SD.set_edge_label(8, 6, 'no label')
801	sage: SD.set_edge_label(13, 17, 'k > h')
802	sage: SD.set_edge_label(13, 14, 'k = h')
803	sage: SD.set_edge_label(17, 15, 'v_k finite')
804	sage: SD.set_edge_label(14, 15, 'v_k m.c.r.')
805	sage: posn = {1:[ 3,-3], 2:[0,2], 3:[0, 13], 4:[3,9], 5:[3,3], 6:[16, 13], 7:[6,1], 8:[6,6], 9:[6,11], 10:[9,1], 11:[10,6], 12:[13,6], 13:[16,2], 14:[10,-6], 15:[0,-10], 16:[14,-6], 17:[16,-10], 18:[6,-4]}
806	sage: SD.plot(pos=posn, vertex_size=400, vertex_colors={'#FFFFFF':range(1,19)}, edge_labels=True).save('search_tree.png')
807
808	EXAMPLES:
809	sage: import sage.graphs.graph_isom
810	sage: from sage.graphs.graph_isom import search_tree
811	sage: from sage.graphs.graph import enum
812	sage: from sage.groups.perm_gps.permgroup import PermutationGroup # long time
813	sage: from sage.graphs.graph_isom import perm_group_elt # long time
814
815	sage: G = graphs.DodecahedralGraph()
816	sage: Pi=[range(20)]
817	sage: a,b = search_tree(G, Pi)
818	sage: print a, enum(b)
819	[[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], [2, 1, 0, 19, 18, 11, 10, 9, 8, 7, 6, 5, 15, 14, 13, 12, 16, 17, 4, 3]] 17318942212009113839976787462421724338461987195898671092180383421848885858584973127639899792828728124797968735273000
820	sage: c = search_tree(G, Pi, lab=False)
821	sage: print c
822	[[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], [2, 1, 0, 19, 18, 11, 10, 9, 8, 7, 6, 5, 15, 14, 13, 12, 16, 17, 4, 3]]
823	sage: DodecAut = PermutationGroup([perm_group_elt(aa) for aa in a]) # long time
824	sage: DodecAut.character_table() # long time
825	[ 1 1 1 1 1 1 1 1 1 1]
826	[ 1 -1 1 1 -1 1 -1 1 -1 -1]
827	[ 3 -1 0 -1 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 3]
828	[ 3 -1 0 -1 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 3]
829	[ 3 1 0 -1 zeta5^3 + zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 -zeta5^3 - zeta5^2 - 1 -3]
830	[ 3 1 0 -1 -zeta5^3 - zeta5^2 - 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 zeta5^3 + zeta5^2 -3]
831	[ 4 0 1 0 -1 -1 1 -1 -1 4]
832	[ 4 0 1 0 1 -1 -1 -1 1 -4]
833	[ 5 1 -1 1 0 0 -1 0 0 5]
834	[ 5 -1 -1 1 0 0 1 0 0 -5]
835	sage: DodecAut2 = PermutationGroup([perm_group_elt(cc) for cc in c]) # long time
836	sage: DodecAut2.character_table() # long time
837	[ 1 1 1 1 1 1 1 1 1 1]
838	[ 1 -1 1 1 -1 1 -1 1 -1 -1]
839	[ 3 -1 0 -1 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 zeta5^3 + zeta5^2 + 1 3]
840	[ 3 -1 0 -1 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 -zeta5^3 - zeta5^2 3]
841	[ 3 1 0 -1 zeta5^3 + zeta5^2 zeta5^3 + zeta5^2 + 1 0 -zeta5^3 - zeta5^2 -zeta5^3 - zeta5^2 - 1 -3]
842	[ 3 1 0 -1 -zeta5^3 - zeta5^2 - 1 -zeta5^3 - zeta5^2 0 zeta5^3 + zeta5^2 + 1 zeta5^3 + zeta5^2 -3]
843	[ 4 0 1 0 -1 -1 1 -1 -1 4]
844	[ 4 0 1 0 1 -1 -1 -1 1 -4]
845	[ 5 1 -1 1 0 0 -1 0 0 5]
846	[ 5 -1 -1 1 0 0 1 0 0 -5]
847
848	sage: G = graphs.PetersenGraph()
849	sage: Pi=[range(10)]
850	sage: a,b = search_tree(G, Pi)
851	sage: print a, enum(b)
852	[[0, 1, 2, 7, 5, 4, 6, 3, 9, 8], [0, 1, 6, 8, 5, 4, 2, 9, 3, 7], [0, 4, 3, 8, 5, 1, 9, 2, 6, 7], [1, 0, 4, 9, 6, 2, 5, 3, 7, 8], [2, 1, 0, 5, 7, 3, 6, 4, 8, 9]] 8715233764864019919698297664
853	sage: c = search_tree(G, Pi, lab=False)
854	sage: PAut = PermutationGroup([perm_group_elt(aa) for aa in a]) # long time
855	sage: PAut.character_table() # long time
856	[ 1 1 1 1 1 1 1]
857	[ 1 -1 1 -1 1 -1 1]
858	[ 4 -2 0 1 1 0 -1]
859	[ 4 2 0 -1 1 0 -1]
860	[ 5 1 1 1 -1 -1 0]
861	[ 5 -1 1 -1 -1 1 0]
862	[ 6 0 -2 0 0 0 1]
863	sage: PAut = PermutationGroup([perm_group_elt(cc) for cc in c]) # long time
864	sage: PAut.character_table() # long time
865	[ 1 1 1 1 1 1 1]
866	[ 1 -1 1 -1 1 -1 1]
867	[ 4 -2 0 1 1 0 -1]
868	[ 4 2 0 -1 1 0 -1]
869	[ 5 1 1 1 -1 -1 0]
870	[ 5 -1 1 -1 -1 1 0]
871	[ 6 0 -2 0 0 0 1]
872
873	sage: G = graphs.CubeGraph(3)
874	sage: Pi = []
875	sage: for i in range(8):
876	... b = Integer(i).binary()
877	... Pi.append(b.zfill(3))
878	...
879	sage: Pi = [Pi]
880	sage: a,b = search_tree(G, Pi)
881	sage: print a, enum(b)
882	[[0, 2, 1, 3, 4, 6, 5, 7], [0, 1, 4, 5, 2, 3, 6, 7], [1, 0, 3, 2, 5, 4, 7, 6]] 520239721777506480
883	sage: c = search_tree(G, Pi, lab=False)
884
885	sage: PermutationGroup([perm_group_elt(aa) for aa in a]).order() # long time
886	48
887	sage: PermutationGroup([perm_group_elt(cc) for cc in c]).order() # long time
888	48
889	sage: DodecAut.order() # long time
890	120
891	sage: PAut.order() # long time
892	120
893
894	sage: D = graphs.DodecahedralGraph()
895	sage: a,b,c = search_tree(D, [range(20)], certify=True)
896	sage: from sage.plot.plot import GraphicsArray # long time
897	sage: import networkx # long time
898	sage: position_D = networkx.spring_layout(D._nxg) # long time
899	sage: position_b = {} # long time
900	sage: for vert in position_D: # long time
901	... position_b[c[vert]] = position_D[vert]
902	sage: graphics_array([D.plot(pos=position_D), b.plot(pos=position_b)]).show() # long time
903	sage: c
904	{0: 0, 1: 19, 2: 16, 3: 15, 4: 9, 5: 1, 6: 10, 7: 8, 8: 14, 9: 12, 10: 17, 11: 11, 12: 5, 13: 6, 14: 2, 15: 4, 16: 3, 17: 7, 18: 13, 19: 18}
905
906	BENCHMARKS:
907	The following examples are given to check modifications to the algorithm
908	for optimization.
909
910	sage: G = Graph({0:[]})
911	sage: Pi = [[0]]
912	sage: a,b = search_tree(G, Pi)
913	sage: print a, enum(b)
914	[] 0
915	sage: a,b = search_tree(G, Pi, dig=True)
916	sage: print a, enum(b)
917	[] 0
918	sage: search_tree(G, Pi, lab=False)
919	[]
920
921	sage: from sage.graphs.graph_isom import all_labeled_graphs, all_ordered_partitions
922
923	sage: graph2 = all_labeled_graphs(2)
924	sage: part2 = all_ordered_partitions(range(2))
925	sage: for G in graph2:
926	... for Pi in part2:
927	... a,b = search_tree(G, Pi)
928	... c,d = search_tree(G, Pi, dig=True)
929	... e = search_tree(G, Pi, lab=False)
930	... a = str(a); b = str(enum(b)); c = str(c); d = str(enum(d)); e = str(e)
931	... print a.ljust(15), b.ljust(5), c.ljust(15), d.ljust(5), e.ljust(15)
932	[] 0 [] 0 []
933	[] 0 [] 0 []
934	[[1, 0]] 0 [[1, 0]] 0 [[1, 0]]
935	[[1, 0]] 0 [[1, 0]] 0 [[1, 0]]
936	[] 6 [] 6 []
937	[] 6 [] 6 []
938	[[1, 0]] 6 [[1, 0]] 6 [[1, 0]]
939	[[1, 0]] 6 [[1, 0]] 6 [[1, 0]]
940
941	sage: graph3 = all_labeled_graphs(3)
942	sage: part3 = all_ordered_partitions(range(3))
943	sage: for G in graph3:
944	... for Pi in part3:
945	... a,b = search_tree(G, Pi)
946	... c,d = search_tree(G, Pi, dig=True)
947	... e = search_tree(G, Pi, lab=False)
948	... a = str(a); b = str(enum(b)); c = str(c); d = str(enum(d)); e = str(e)
949	... print a.ljust(15), b.ljust(5), c.ljust(15), d.ljust(5), e.ljust(15)
950	[] 0 [] 0 []
951	[] 0 [] 0 []
952	[[0, 2, 1]] 0 [[0, 2, 1]] 0 [[0, 2, 1]]
953	[[0, 2, 1]] 0 [[0, 2, 1]] 0 [[0, 2, 1]]
954	[] 0 [] 0 []
955	[] 0 [] 0 []
956	[[2, 1, 0]] 0 [[2, 1, 0]] 0 [[2, 1, 0]]
957	[[2, 1, 0]] 0 [[2, 1, 0]] 0 [[2, 1, 0]]
958	[] 0 [] 0 []
959	[] 0 [] 0 []
960	[[1, 0, 2]] 0 [[1, 0, 2]] 0 [[1, 0, 2]]
961	[[1, 0, 2]] 0 [[1, 0, 2]] 0 [[1, 0, 2]]
962	[[1, 0, 2]] 0 [[1, 0, 2]] 0 [[1, 0, 2]]
963	[[2, 1, 0]] 0 [[2, 1, 0]] 0 [[2, 1, 0]]
964	[[1, 0, 2]] 0 [[1, 0, 2]] 0 [[1, 0, 2]]
965	[[0, 2, 1]] 0 [[0, 2, 1]] 0 [[0, 2, 1]]
966	[[2, 1, 0]] 0 [[2, 1, 0]] 0 [[2, 1, 0]]
967	[[0, 2, 1]] 0 [[0, 2, 1]] 0 [[0, 2, 1]]
968	[[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]]
969	[[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]]
970	[[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]]
971	[[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]]
972	[[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]]
973	[[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]] 0 [[0, 2, 1], [1, 0, 2]]
974	[] 10 [] 10 []
975	[] 10 [] 10 []
976	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
977	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
978	[] 68 [] 68 []
979	[] 160 [] 160 []
980	[] 68 [] 68 []
981	[] 68 [] 68 []
982	[] 68 [] 68 []
983	[] 160 [] 160 []
984	[] 68 [] 68 []
985	[] 68 [] 68 []
986	[] 10 [] 10 []
987	[] 10 [] 10 []
988	[] 10 [] 10 []
989	[[0, 2, 1]] 160 [[0, 2, 1]] 160 [[0, 2, 1]]
990	[] 10 [] 10 []
991	[[0, 2, 1]] 160 [[0, 2, 1]] 160 [[0, 2, 1]]
992	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
993	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
994	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
995	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
996	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
997	[[0, 2, 1]] 10 [[0, 2, 1]] 10 [[0, 2, 1]]
998	[] 68 [] 68 []
999	[] 160 [] 160 []
1000	[] 68 [] 68 []
1001	[] 68 [] 68 []
1002	[] 10 [] 10 []
1003	[] 10 [] 10 []
1004	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1005	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1006	[] 160 [] 160 []
1007	[] 68 [] 68 []
1008	[] 68 [] 68 []
1009	[] 68 [] 68 []
1010	[] 10 [] 10 []
1011	[[2, 1, 0]] 160 [[2, 1, 0]] 160 [[2, 1, 0]]
1012	[] 10 [] 10 []
1013	[] 10 [] 10 []
1014	[[2, 1, 0]] 160 [[2, 1, 0]] 160 [[2, 1, 0]]
1015	[] 10 [] 10 []
1016	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1017	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1018	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1019	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1020	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1021	[[2, 1, 0]] 10 [[2, 1, 0]] 10 [[2, 1, 0]]
1022	[] 78 [] 78 []
1023	[] 170 [] 170 []
1024	[] 78 [] 78 []
1025	[] 78 [] 78 []
1026	[] 78 [] 78 []
1027	[] 170 [] 170 []
1028	[] 78 [] 78 []
1029	[] 78 [] 78 []
1030	[] 228 [] 228 []
1031	[] 228 [] 228 []
1032	[[1, 0, 2]] 228 [[1, 0, 2]] 228 [[1, 0, 2]]
1033	[[1, 0, 2]] 228 [[1, 0, 2]] 228 [[1, 0, 2]]
1034	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1035	[] 170 [] 170 []
1036	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1037	[] 170 [] 170 []
1038	[] 170 [] 170 []
1039	[] 170 [] 170 []
1040	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1041	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1042	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1043	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1044	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1045	[[1, 0, 2]] 78 [[1, 0, 2]] 78 [[1, 0, 2]]
1046	[] 160 [] 160 []
1047	[] 68 [] 68 []
1048	[] 68 [] 68 []
1049	[] 68 [] 68 []
1050	[] 160 [] 160 []
1051	[] 68 [] 68 []
1052	[] 68 [] 68 []
1053	[] 68 [] 68 []
1054	[] 10 [] 10 []
1055	[] 10 [] 10 []
1056	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1057	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1058	[[1, 0, 2]] 160 [[1, 0, 2]] 160 [[1, 0, 2]]
1059	[] 10 [] 10 []
1060	[[1, 0, 2]] 160 [[1, 0, 2]] 160 [[1, 0, 2]]
1061	[] 10 [] 10 []
1062	[] 10 [] 10 []
1063	[] 10 [] 10 []
1064	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1065	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1066	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1067	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1068	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1069	[[1, 0, 2]] 10 [[1, 0, 2]] 10 [[1, 0, 2]]
1070	[] 170 [] 170 []
1071	[] 78 [] 78 []
1072	[] 78 [] 78 []
1073	[] 78 [] 78 []
1074	[] 228 [] 228 []
1075	[] 228 [] 228 []
1076	[[2, 1, 0]] 228 [[2, 1, 0]] 228 [[2, 1, 0]]
1077	[[2, 1, 0]] 228 [[2, 1, 0]] 228 [[2, 1, 0]]
1078	[] 78 [] 78 []
1079	[] 170 [] 170 []
1080	[] 78 [] 78 []
1081	[] 78 [] 78 []
1082	[] 170 [] 170 []
1083	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1084	[] 170 [] 170 []
1085	[] 170 [] 170 []
1086	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1087	[] 170 [] 170 []
1088	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1089	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1090	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1091	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1092	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1093	[[2, 1, 0]] 78 [[2, 1, 0]] 78 [[2, 1, 0]]
1094	[] 228 [] 228 []
1095	[] 228 [] 228 []
1096	[[0, 2, 1]] 228 [[0, 2, 1]] 228 [[0, 2, 1]]
1097	[[0, 2, 1]] 228 [[0, 2, 1]] 228 [[0, 2, 1]]
1098	[] 170 [] 170 []
1099	[] 78 [] 78 []
1100	[] 78 [] 78 []
1101	[] 78 [] 78 []
1102	[] 170 [] 170 []
1103	[] 78 [] 78 []
1104	[] 78 [] 78 []
1105	[] 78 [] 78 []
1106	[] 170 [] 170 []
1107	[] 170 [] 170 []
1108	[] 170 [] 170 []
1109	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1110	[] 170 [] 170 []
1111	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1112	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1113	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1114	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1115	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1116	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1117	[[0, 2, 1]] 78 [[0, 2, 1]] 78 [[0, 2, 1]]
1118	[] 238 [] 238 []
1119	[] 238 [] 238 []
1120	[[0, 2, 1]] 238 [[0, 2, 1]] 238 [[0, 2, 1]]
1121	[[0, 2, 1]] 238 [[0, 2, 1]] 238 [[0, 2, 1]]
1122	[] 238 [] 238 []
1123	[] 238 [] 238 []
1124	[[2, 1, 0]] 238 [[2, 1, 0]] 238 [[2, 1, 0]]
1125	[[2, 1, 0]] 238 [[2, 1, 0]] 238 [[2, 1, 0]]
1126	[] 238 [] 238 []
1127	[] 238 [] 238 []
1128	[[1, 0, 2]] 238 [[1, 0, 2]] 238 [[1, 0, 2]]
1129	[[1, 0, 2]] 238 [[1, 0, 2]] 238 [[1, 0, 2]]
1130	[[1, 0, 2]] 238 [[1, 0, 2]] 238 [[1, 0, 2]]
1131	[[2, 1, 0]] 238 [[2, 1, 0]] 238 [[2, 1, 0]]
1132	[[1, 0, 2]] 238 [[1, 0, 2]] 238 [[1, 0, 2]]
1133	[[0, 2, 1]] 238 [[0, 2, 1]] 238 [[0, 2, 1]]
1134	[[2, 1, 0]] 238 [[2, 1, 0]] 238 [[2, 1, 0]]
1135	[[0, 2, 1]] 238 [[0, 2, 1]] 238 [[0, 2, 1]]
1136	[[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]]
1137	[[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]]
1138	[[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]]
1139	[[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]]
1140	[[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]]
1141	[[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]] 238 [[0, 2, 1], [1, 0, 2]]
1142
1143	sage: C = graphs.CubeGraph(1)
1144	sage: gens = search_tree(C, [C.vertices()], lab=False)
1145	sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() # long time
1146	2
1147	sage: C = graphs.CubeGraph(2)
1148	sage: gens = search_tree(C, [C.vertices()], lab=False)
1149	sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() # long time
1150	8
1151	sage: C = graphs.CubeGraph(3)
1152	sage: gens = search_tree(C, [C.vertices()], lab=False)
1153	sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() # long time
1154	48
1155	sage: C = graphs.CubeGraph(4)
1156	sage: gens = search_tree(C, [C.vertices()], lab=False)
1157	sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() # long time
1158	384
1159	sage: C = graphs.CubeGraph(5)
1160	sage: gens = search_tree(C, [C.vertices()], lab=False)
1161	sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() # long time
1162	3840
1163	sage: C = graphs.CubeGraph(6)
1164	sage: gens = search_tree(C, [C.vertices()], lab=False)
1165	sage: PermutationGroup([perm_group_elt(aa) for aa in gens]).order() # long time
1166	46080
1167	"""
1168	cdef int i, j, # local variables
1169
1170	cdef OrbitPartition Theta, OP
1171	cdef int index = 0, size = 1 # see Theorem 2.33 in [1]
1172
1173	cdef int L = 100 # memory limit for storing values from fix and mcr:
1174	# Phi and Omega store specific information about some
1175	# of the automorphisms we already know about, and they
1176	# are arrays of length L
1177	cdef int **Phi # stores the fixed point sets of each automorphism
1178	cdef int **Omega # stores the minimal elements of each cell of the
1179	# orbit partition
1180	cdef int l = -1 # current index for storing values in Phi and Omega-
1181	# we start at -1 so that when we increment first,
1182	# the first place we write to is 0.
1183	cdef int **W # for each k, W[k] is a list of the vertices to be searched
1184	# down from the current partition nest, at k
1185	# Phi and Omega are ultimately used to make the size of W
1186	# as small as possible
1187
1188	cdef PartitionStack nu, zeta, rho
1189	cdef int k_rho # the number of partitions in rho
1190	cdef int k = 0 # the number of partitions in nu
1191	cdef int h = -1 # longest common ancestor of zeta and nu:
1192	# zeta[h] == nu[h], zeta[h+1] != nu[h+1]
1193	cdef int hb # longest common ancestor of rho and nu:
1194	# rho[hb] == nu[hb], rho[hb+1] != nu[hb+1]
1195	cdef int hh = 1 # the height of the oldest ancestor of nu
1196	# satisfying Lemma 2.25 in [1]
1197	cdef int ht # smallest such that all descendants of zeta[ht] are equivalent
1198
1199	cdef mpz_t Lambda_mpz, zf_mpz, *zb_mpz # for tracking indicator values
1200	# zf and zb are indicator vectors remembering Lambda[k] for zeta and rho,
1201	# respectively
1202	cdef int hzf # the max height for which Lambda and zf agree
1203	cdef int hzb = -1 # the max height for which Lambda and zb agree
1204
1205	cdef int **M # for the square adjacency matrix
1206	cdef int *gamma # for storing permutations
1207	cdef int *alpha # for storing pointers to cells of nu[k]:
1208	# allocated to be length 4*n for scratch (see functions
1209	# _sort_by_function and _refine_by_square_matrix)
1210	cdef int *v # list of vertices determining nu
1211	cdef int *e # 0 or 1, see states 12 and 17
1212	cdef int state # keeps track of place in algorithm
1213	cdef int _dig, tvc, tvh, n = G.order()
1214
1215	# trivial case
1216	if n == 0:
1217	if lab:
1218	H = G.copy()
1219	if dict:
1220	ddd = {}
1221	if certify:
1222	dd = {}
1223	if dict:
1224	return [[]], ddd, H, dd
1225	else:
1226	return [[]], H, dd
1227	if lab and dict:
1228	return [[]], ddd, H
1229	elif lab:
1230	return [[]], H
1231	elif dict:
1232	return [[]], ddd
1233	else:
1234	return [[]]
1235
1236	# create to and from mappings to relabel vertices to the set {0,...,n-1}
1237	listto = G.vertices()
1238	ffrom = {}
1239	for vvv in listto:
1240	ffrom[vvv] = listto.index(vvv)
1241	to = {}
1242	for i from 0 <= i < len(listto):
1243	to[i] = listto[i]
1244	G.relabel(ffrom)
1245	Pi2 = []
1246	for cell in Pi:
1247	Pi2.append([ffrom[c] for c in cell])
1248	Pi = Pi2
1249
1250	# allocate pointers
1251	W = <int *> sage_malloc( 2 (n + L) * sizeof(int *) )
1252	if not W:
1253	raise MemoryError("Error allocating memory. Perhaps you are out?")
1254	M = W + n
1255	Phi = W + 2*n
1256	Omega = W + 2*n + L
1257
1258	# allocate pointers for GMP ints
1259	Lambda_mpz = <mpz_t > sage_malloc( 3 (n+2) * sizeof(mpz_t) )
1260	if not Lambda_mpz:
1261	sage_free(W)
1262	raise MemoryError("Error allocating memory. Perhaps you are out?")
1263	zf_mpz = Lambda_mpz + n + 2
1264	zb_mpz = Lambda_mpz + 2*n + 4
1265
1266	# allocate arrays
1267	gamma = <int > sage_malloc( ( n ( 2 * (n + L) + 7 ) + 1 ) * sizeof(int) )
1268	if not gamma:
1269	sage_free(W)
1270	sage_free(Lambda_mpz)
1271	raise MemoryError("Error allocating memory. Perhaps you are out?")
1272	alpha = gamma + n( 2(n + L) + 1 )
1273	v = alpha + 4*n + 1
1274	e = v + n
1275	for i from 0 <= i < n:
1276	W[i] = gamma + n + n*i
1277	M[i] = gamma + n( 1 + n + 2L + i )
1278
1279	# allocate GMP ints
1280	for i from 0 <= i < n+2:
1281	mpz_init(Lambda_mpz[i])
1282	mpz_init_set_si(zf_mpz[i], -1) # correspond to default values of
1283	mpz_init_set_si(zb_mpz[i], -1) # "infinity"
1284	# Note that there is a potential memory leak here - if a particular
1285	# mpz fails to allocate, this is not checked for
1286	for i from 0 <= i < L:
1287	Phi[i] = gamma + n*( 1 + n + i )
1288	Omega[i] = gamma + n*( 1 + n + i + L )
1289
1290	# create the dense boolean matrix
1291	for i from 0 <= i < n:
1292	for j from 0 <= j < n:
1293	M[i][j] = 0
1294	W[i][j] = 0
1295	if isinstance(G, Graph):
1296	for i, j, la in G.edge_iterator():
1297	M[i][j] = 1
1298	M[j][i] = 1
1299	elif isinstance(G, DiGraph):
1300	for i, j, la in G.edge_iterator():
1301	M[i][j] = 1
1302
1303	# set up the rest of the variables
1304	nu = PartitionStack(Pi)
1305	Theta = OrbitPartition(n)
1306	G_enum = _enumerate_graph(M, n)
1307	output = []
1308	if dig: _dig = 1
1309	else: _dig = 0
1310	if certify:
1311	lab=True
1312
1313	if verbosity > 1:
1314	t = cputime()
1315	if verbosity > 2:
1316	rho = PartitionStack(n)
1317	zeta = PartitionStack(n)
1318	state = 1
1319	while state != -1:
1320	if verbosity > 0:
1321	print '-----'
1322	print 'state:', state
1323	print 'nu'
1324	print [nu.entries[iii] for iii in range(n)]
1325	print [nu.levels[iii] for iii in range(n)]
1326	if verbosity > 1:
1327	t = cputime(t)
1328	print 'time:', t
1329	if verbosity > 2:
1330	print 'k: ' + str(k)
1331	print 'zeta:'
1332	print [zeta.entries[iii] for iii in range(n)]
1333	print [zeta.levels[iii] for iii in range(n)]
1334	print 'rho'
1335	print [rho.entries[iii] for iii in range(n)]
1336	print [rho.levels[iii] for iii in range(n)]
1337	if verbosity > 3:
1338	Thetarep = []
1339	for i from 0 <= i < n:
1340	j = Theta._find(i)
1341	didit = False
1342	for celll in Thetarep:
1343	if celll[0] == j:
1344	celll.append(i)
1345	didit = True
1346	if not didit:
1347	Thetarep.append([j])
1348	print 'Theta: ', str(Thetarep)
1349
1350	if state == 1: # Entry point to algorithm
1351	# get alpha to point to cells of nu
1352	j = 1
1353	alpha[0] = 0
1354	for i from 0 < i < n:
1355	if nu.levels[i-1] == 0:
1356	alpha[j] = i
1357	j += 1
1358	alpha[j] = -1
1359
1360	# "nu[0] := R(G, Pi, Pi)"
1361	nu._refine_by_square_matrix(k, alpha, n, M, _dig)
1362
1363	if not _dig:
1364	if nu._sat_225(k, n): hh = k
1365	if nu._is_discrete(k): state = 18; continue
1366
1367	# store the first smallest nontrivial cell in W[k], and set v[k]
1368	# equal to its minimum element
1369	v[k] = nu._first_smallest_nontrivial(W[k], k, n)
1370	mpz_set_ui(Lambda_mpz[k], 0)
1371	e[k] = 0 # see state 12, and 17
1372	state = 2
1373
1374	elif state == 2: # Move down the search tree one level by refining nu
1375	k += 1
1376
1377	# "nu[k] := nu[k-1] perp v[k-1]"
1378	nu._clear(k)
1379	alpha[0] = nu._split_vertex(v[k-1], k)
1380	alpha[1] = -1
1381	i = nu._refine_by_square_matrix(k, alpha, n, M, _dig)
1382
1383	# add one, then multiply by the invariant
1384	mpz_add_ui(Lambda_mpz[k], Lambda_mpz[k-1], 1)
1385	mpz_mul_si(Lambda_mpz[k], Lambda_mpz[k], i)
1386
1387	# only if this is the first time moving down the search tree:
1388	if h == -1: state = 5; continue
1389
1390	# update hzf
1391	if hzf == k-1 and mpz_cmp(Lambda_mpz[k], zf_mpz[k]) == 0: hzf = k
1392	if not lab: state = 3; continue
1393
1394	# "qzb := cmp(Lambda[k], zb[k])"
1395	if mpz_cmp_si(zb_mpz[k], -1) == 0: # if "zb[k] == oo"
1396	qzb = -1
1397	else:
1398	qzb = mpz_cmp( Lambda_mpz[k], zb_mpz[k] )
1399	# update hzb
1400	if hzb == k-1 and qzb == 0: hzb = k
1401
1402	# if Lambda[k] > zb[k], then zb[k] := Lambda[k]
1403	# (zb keeps track of the indicator invariants corresponding to
1404	# rho, the closest canonical leaf so far seen- if Lambda is
1405	# bigger, then rho must be about to change
1406	if qzb > 0: mpz_set(zb_mpz[k], Lambda_mpz[k])
1407	state = 3
1408
1409	elif state == 3: # attempt to rule out automorphisms while moving down
1410	# the tree
1411	if hzf <= k or (lab and qzb >= 0): # changed hzb to hzf, == to <=
1412	state = 4
1413	else: state = 6
1414	# if k > hzf, then we know that nu currently does not look like
1415	# zeta, the first terminal node encountered. Then if we are not
1416	# looking for a canonical label, there is no reason to continue.
1417	# However, if we are looking for one, and qzb < 0, i.e.
1418	# Lambda[k] < zb[k], then the indicator is not maximal, and we
1419	# can't reach a canonical leaf.
1420
1421	elif state == 4: # at this point we have -not- ruled out the presence
1422	# of automorphisms
1423	if nu._is_discrete(k): state = 7; continue
1424
1425	# store the first smallest nontrivial cell in W[k], and set v[k]
1426	# equal to its minimum element
1427	v[k] = nu._first_smallest_nontrivial(W[k], k, n)
1428
1429	if _dig or not nu._sat_225(k, n): hh = k + 1
1430	e[k] = 0 # see state 12, and 17
1431	state = 2 # continue down the tree
1432
1433	elif state == 5: # alternative to 3: since we have not yet gotten
1434	# zeta, there are no automorphisms to rule out.
1435	# instead we record Lambda to zf and zb
1436	# (see state 3)
1437	mpz_set(zf_mpz[k], Lambda_mpz[k])
1438	mpz_set(zb_mpz[k], Lambda_mpz[k])
1439	state = 4
1440
1441	elif state == 6: # at this stage, there is no reason to continue
1442	# downward, and an automorphism has not been
1443	# discovered
1444	j = k
1445
1446	# return to the longest ancestor nu[i] of nu that could have a
1447	# descendant equivalent to zeta or could improve on rho.
1448	# All terminal nodes descending from nu[hh] are known to be
1449	# equivalent, so i < hh. Also, if i > hzb, none of the
1450	# descendants of nu[i] can improve rho, since the indicator is
1451	# off (Lambda(nu) < Lambda(rho)). If i >= ht, then no descendant
1452	# of nu[i] is equivalent to zeta (see [1, p67]).
1453	if ht-1 > hzb:
1454	if ht-1 < hh-1:
1455	k = ht-1
1456	else:
1457	k = hh-1
1458	else:
1459	if hzb < hh-1:
1460	k = hzb
1461	else:
1462	k = hh-1
1463
1464	# TODO: investigate the following line
1465	if k == -1: k = 0 # not in BDM, broke at G = Graph({0:[], 1:[]}), Pi = [[0,1]], lab=False
1466
1467	if j == hh: state = 13; continue
1468	# recall hh: the height of the oldest ancestor of zeta for which
1469	# Lemma 2.25 is satsified, which implies that all terminal nodes
1470	# descended from there are equivalent (or simply k if 2.25 does
1471	# not apply). If we are looking at such a node, then the partition
1472	# at nu[hh] can be used for later pruning, so we store its fixed
1473	# set and a set of representatives of its cells
1474	if l < L-1: l += 1
1475	for i from 0 <= i < n:
1476	Omega[l][i] = 0 # changed Lambda to Omega
1477	Phi[l][i] = 0
1478	if nu._is_min_cell_rep(i, hh):
1479	Omega[l][i] = 1
1480	if nu._is_fixed(i, hh):
1481	Phi[l][i] = 1
1482
1483	state = 12
1484
1485	elif state == 7: # we have just arrived at a terminal node of the
1486	# search tree T(G, Pi)
1487	# if this is the first terminal node, go directly to 18, to
1488	# process zeta
1489	if h == -1: state = 18; continue
1490
1491	# hzf is the extremal height of ancestors of both nu and zeta,
1492	# so if k < hzf, nu is not equivalent to zeta, i.e. there is no
1493	# automorphism to discover.
1494	# TODO: investigate why, in practice, the same does not seem to be
1495	# true for hzf < k... BDM had !=, not <, and this broke at
1496	# G = Graph({0:[],1:[],2:[]}), Pi = [[0,1,2]]
1497	if k < hzf: state = 8; continue
1498
1499	# get the permutation corresponding to this terminal node
1500	nu._get_permutation_from(zeta, gamma)
1501
1502	if verbosity > 3:
1503	print 'automorphism discovered:'
1504	print [gamma[iii] for iii in range(n)]
1505
1506	# if G^gamma == G, the permutation is an automorphism, goto 10
1507	if G_enum == _enumerate_graph_with_permutation(M, n, gamma):
1508	state = 10
1509	else:
1510	state = 8
1511
1512	elif state == 8: # we have just ruled out the presence of automorphism
1513	# and have not yet considered whether nu improves on
1514	# rho
1515	# if we are not searching for a canonical label, there is nothing
1516	# to do here
1517	if (not lab) or (qzb < 0):
1518	state = 6; continue
1519
1520	# if Lambda[k] > zb[k] or nu is shorter than rho, then we have
1521	# found an improvement for rho
1522	if (qzb > 0) or (k < k_rho):
1523	state = 9; continue
1524
1525	# if G(nu) > G(rho), goto 9
1526	# if G(nu) < G(rho), goto 6
1527	# if G(nu) == G(rho), get the automorphism and goto 10
1528	m1 = nu._enumerate_graph_from_discrete(M, n)
1529	m2 = rho._enumerate_graph_from_discrete(M, n)
1530
1531	if m1 > m2:
1532	state = 9; continue
1533	if m1 < m2:
1534	state = 6; continue
1535
1536	rho._get_permutation_from(nu, gamma)
1537	if verbosity > 3:
1538	print 'automorphism discovered:'
1539	print [gamma[iii] for iii in range(n)]
1540	state = 10
1541
1542	elif state == 9: # entering this state, nu is a best-so-far guess at
1543	# the canonical label
1544	rho = PartitionStack(nu)
1545	k_rho = k
1546
1547	qzb = 0
1548	hb = k
1549	hzb = k
1550
1551	# set zb[k+1] = Infinity
1552	mpz_set_si(zb_mpz[k+1], -1)
1553	state = 6
1554
1555	elif state == 10: # we have an automorphism to process
1556	# increment l
1557	if l < L - 1:
1558	l += 1
1559
1560	# retrieve the orbit partition, and record the relevant
1561	# information
1562	# TODO: this step could be optimized. The variable OP is not
1563	# really necessary
1564	OP = _orbit_partition_from_list_perm(gamma, n)
1565	for i from 0 <= i < n:
1566	Omega[l][i] = OP._is_min_cell_rep(i)
1567	Phi[l][i] = OP._is_fixed(i)
1568
1569	# if each orbit of gamma is part of an orbit in Theta, then the
1570	# automorphism is already in the span of those we have seen
1571	if OP._is_finer_than(Theta, n):
1572	state = 11
1573	continue
1574	# otherwise, incorporate this into Theta
1575	Theta._vee_with(OP, n)
1576
1577	# record the automorphism
1578	output.append([ Integer(gamma[i]) for i from 0 <= i < n ])
1579
1580	# The variable tvc was set to be the minimum element of W[k]
1581	# the last time we were at state 13 and at a node descending to
1582	# zeta. If this is a minimal cell representative of Theta and
1583	# we are searching for a canonical label, goto state 11, i.e.
1584	# backtrack to the common ancestor of rho and nu, then goto state
1585	# 12, i.e. consider whether we still need to search downward from
1586	# there. TODO: explain why
1587	if Theta.elements[tvc] == -1 and lab: ## added "and lab"
1588	state = 11
1589	continue
1590	k = h
1591	state = 13
1592
1593	elif state == 11: # if we are searching for a label, backtrack to the
1594	# common ancestor of nu and rho
1595	k = hb
1596	state = 12
1597
1598	elif state == 12: # we are looking at a branch we may have to continue
1599	# to search downward on
1600	# e keeps track of the motion through the search tree. It is set to
1601	# 1 when you have just finished coming up the search tree, and are
1602	# at a node in the tree for which there may be more branches left
1603	# to explore. In this case, intersect W[k] with Omega[l], since
1604	# there may be an automorphism mapping one element of W[k] to
1605	# another, hence only one must be investigated.
1606	if e[k] == 1:
1607	for j from 0 <= j < n:
1608	if W[k][j] and not Omega[l][j]:
1609	W[k][j] = 0
1610	state = 13
1611
1612	elif state == 13: # hub state
1613	if k == -1:
1614	# the algorithm has finished
1615	state = -1; continue
1616	if k > h:
1617	# if we are not at a node of zeta
1618	state = 17; continue
1619	if k == h:
1620	# if we are at a node of zeta, then state 14 can rule out
1621	# vertices to consider
1622	state = 14; continue
1623
1624	# thus, it must be that k < h, and this means we are done
1625	# searching underneath zeta[k+1], so now, k is the new longest
1626	# ancestor of nu and zeta:
1627	h = k
1628
1629	# set tvc and tvh to the minimum cell representative of W[k]
1630	# (see states 10 and 14)
1631	for i from 0 <= i < n:
1632	if W[k][i]:
1633	tvc = i
1634	break
1635	tvh = tvc
1636	state = 14
1637
1638	elif state == 14: # iterate v[k] through W[k] until a minimum cell rep
1639	# of Theta is found
1640	# The variable tvh was set to be the minimum element of W[k]
1641	# the last time we were at state 13 and at a node descending to
1642	# zeta. If this is in the same cell of Theta as v[k], increment
1643	# index (see Theorem 2.33 in [1])
1644	if Theta._find(v[k]) == Theta._find(tvh):
1645	index += 1
1646
1647	# find the next v[k] in W[k]
1648	i = v[k] + 1
1649	while i < n and not W[k][i]:
1650	i += 1
1651	if i < n:
1652	v[k] = i
1653	else:
1654	# there is no new vertex to consider at this level
1655	v[k] = -1
1656	state = 16
1657	continue
1658
1659	# if the new v[k] is not a minimum cell representative of Theta,
1660	# then we already considered that rep., and that subtree was
1661	# isomorphic to the one corresponding to v[k]
1662	if Theta.elements[v[k]] != -1: state = 14
1663	else:
1664	# otherwise, we do have a vertex to consider
1665	state = 15
1666
1667	elif state == 15: # we have a new vertex, v[k], that we must split on
1668	# hh is smallest such that nu[hh] satisfies Lemma 2.25. If it is
1669	# larger than k+1, it must be modified, since we are changing that
1670	# part
1671	if k + 1 < hh:
1672	hh = k + 1
1673	# hzf is maximal such that indicators line up for nu and zeta
1674	if k < hzf:
1675	hzf = k
1676	if not lab or hb < k: # changed hzb to hb
1677	# in either case there is no need to update hb, which is the
1678	# length of the common ancestor of nu and rho
1679	state = 2; continue
1680	hb = k # changed hzb to hb
1681	qzb = 0
1682	state = 2
1683
1684	elif state == 16: # backtrack one level in the search tree, recording
1685	# information relevant to Theorem 2.33
1686	j = 0
1687	for i from 0 <= i < n:
1688	if W[k][i]: j += 1
1689	if j == index and ht == k+1: ht = k
1690	size = size*index
1691	index = 0
1692	k -= 1
1693	state = 13
1694
1695	elif state == 17: # you have just finished coming up the search tree,
1696	# and must now consider going back down.
1697	if e[k] == 0:
1698	# intersect W[k] with each Omega[i] such that {v_0,...,v_(k-1)}
1699	# is contained in Phi[i]
1700	for i from 0 <= i <= l:
1701	# check if {v_0,...,v_(k-1)} is contained in Phi[i]
1702	# i.e. fixed pointwise by the automorphisms so far seen
1703	j = 0
1704	while j < k and Phi[i][v[j]]:
1705	j += 1
1706	# if so, only check the minimal orbit representatives
1707	if j == k:
1708	for j from 0 <= j < n:
1709	if W[k][j] and not Omega[i][j]:
1710	W[k][j] = 0
1711	e[k] = 1 # see state 12
1712
1713	# see if there is a relevant vertex to split on:
1714	i = v[k] + 1
1715	while i < n and not W[k][i]:
1716	i += 1
1717	if i < n:
1718	v[k] = i
1719	state = 15
1720	continue
1721	else:
1722	v[k] = -1
1723
1724	# otherwise backtrack one level
1725	k -= 1
1726	state = 13
1727
1728	elif state == 18: # The first time we encounter a terminal node, we
1729	# come straight here to set up zeta. This is a one-
1730	# time state.
1731	# initialize counters for zeta:
1732	h = k # zeta[h] == nu[h]
1733	ht = k # nodes descended from zeta[ht] are all equivalent
1734	hzf = k # max such that indicators for zeta and nu agree
1735
1736	zeta = PartitionStack(nu)
1737
1738	k -= 1
1739	if not lab: state = 13; continue
1740
1741	rho = PartitionStack(nu)
1742
1743	# initialize counters for rho:
1744	k_rho = k # number of partitions in rho
1745	hzb = k # max such that indicators for rho and nu agree - BDM had k+1
1746	hb = k # rho[hb] == nu[hb] - BDM had k+1
1747
1748	qzb = 0 # Lambda[k] == zb[k], so...
1749	state = 13
1750
1751	# deallocate the MP integers
1752	for i from 0 <= i < n:
1753	mpz_clear(Lambda_mpz[i])
1754	mpz_clear(zf_mpz[i])
1755	mpz_clear(zb_mpz[i])
1756
1757	sage_free(W)
1758	sage_free(Lambda_mpz)
1759	sage_free(gamma)
1760
1761	# use to and from mappings to relabel vertices back from the set {0,...,n-1}
1762	if lab:
1763	H = _term_pnest_graph(G, rho)
1764	G.relabel(to)
1765	if dict:
1766	ddd = {}
1767	for vvv in G.vertices(): # v is a C variable
1768	if ffrom[vvv] != 0:
1769	ddd[vvv] = ffrom[vvv]
1770	else:
1771	ddd[vvv] = n
1772
1773	# prepare output
1774	if certify:
1775	dd = {}
1776	for i from 0 <= i < n:
1777	dd[rho.entries[i]] = i
1778	# NOTE - this should take the relabeling into account!
1779	if dict:
1780	return output, ddd, H, dd
1781	else:
1782	return output, H, dd
1783	if lab and dict:
1784	return output, ddd, H
1785	elif lab:
1786	return output, H
1787	elif dict:
1788	return output, ddd
1789	else:
1790	return output
1791
1792	# Benchmarking functions
1793
1794	def all_labeled_graphs(n):
1795	"""
1796	Returns all labeled graphs on n vertices {0,1,...,n-1}. Used in
1797	classifying isomorphism types (naive approach), and more importantly
1798	in benchmarking the search algorithm.
1799
1800	EXAMPLE:
1801	sage: import sage.graphs.graph_isom
1802	sage: from sage.graphs.graph_isom import search_tree, all_labeled_graphs
1803	sage: from sage.graphs.graph import enum
1804	sage: Glist = {}
1805	sage: Giso = {}
1806	sage: for n in range(1,5):
1807	... Glist[n] = all_labeled_graphs(n)
1808	... Giso[n] = []
1809	... for g in Glist[n]:
1810	... a, b = search_tree(g, [range(n)])
1811	... inn = False
1812	... for gi in Giso[n]:
1813	... if enum(b) == enum(gi):
1814	... inn = True
1815	... if not inn:
1816	... Giso[n].append(b)
1817	sage: for n in Giso:
1818	... print n, len(Giso[n])
1819	1 1
1820	2 2
1821	3 4
1822	4 11
1823	sage: n = 5 # long time
1824	sage: Glist[n] = all_labeled_graphs(n) # long time
1825	sage: Giso[n] = [] # long time
1826	sage: for g in Glist[5]: # long time
1827	... a, b = search_tree(g, [range(n)])
1828	... inn = False
1829	... for gi in Giso[n]:
1830	... if enum(b) == enum(gi):
1831	... inn = True
1832	... if not inn:
1833	... Giso[n].append(b)
1834	sage: print n, len(Giso[n]) # long time
1835	5 34
1836	sage: graphs_list.show_graphs(Giso[4])
1837	"""
1838	TE = []
1839	for i in range(n):
1840	for j in range(i):
1841	TE.append((i, j))
1842	m = len(TE)
1843	Glist= []
1844	for i in range(2**m):
1845	G = Graph()
1846	G.add_vertices(range(n))
1847	b = Integer(i).binary()
1848	b = '0'*(m-len(b)) + b
1849	for i in range(m):
1850	if int(b[i]):
1851	G.add_edge(TE[i])
1852	Glist.append(G)
1853	return Glist
1854
1855	def kpow(listy, k):
1856	"""
1857	Returns the subset of the power set of listy consisting of subsets of size
1858	k. Used in all_ordered_partitions.
1859	"""
1860	list = []
1861	if k > 1:
1862	for LL in kpow(listy, k-1):
1863	for a in listy:
1864	if not a in LL:
1865	list.append([a] + LL)
1866	if k == 1:
1867	for i in listy:
1868	list.append([i])
1869	return list
1870
1871	def all_ordered_partitions(listy):
1872	"""
1873	Returns all ordered partitions of the set {0,1,...,n-1}. Used in
1874	benchmarking the search algorithm.
1875	"""
1876	LL = []
1877	for i in range(1,len(listy)+1):
1878	for cell in kpow(listy, i):
1879	list_remainder = [x for x in listy if x not in cell]
1880	remainder_partitions = all_ordered_partitions(list_remainder)
1881	for remainder in remainder_partitions:
1882	LL.append( [cell] + remainder )
1883	if len(listy) == 0:
1884	return [[]]
1885	else:
1886	return LL
1887
1888	def all_labeled_digraphs_with_loops(n):
1889	"""
1890	Returns all labeled digraphs on n vertices {0,1,...,n-1}. Used in
1891	classifying isomorphism types (naive approach), and more importantly
1892	in benchmarking the search algorithm.
1893
1894	EXAMPLE:
1895	sage: import sage.graphs.graph_isom
1896	sage: from sage.graphs.graph_isom import search_tree, all_labeled_digraphs_with_loops
1897	sage: from sage.graphs.graph import enum
1898	sage: Glist = {}
1899	sage: Giso = {}
1900	sage: for n in range(1,4):
1901	... Glist[n] = all_labeled_digraphs_with_loops(n)
1902	... Giso[n] = []
1903	... for g in Glist[n]:
1904	... a, b = search_tree(g, [range(n)], dig=True)
1905	... inn = False
1906	... for gi in Giso[n]:
1907	... if enum(b) == enum(gi):
1908	... inn = True
1909	... if not inn:
1910	... Giso[n].append(b)
1911	sage: for n in Giso:
1912	... print n, len(Giso[n])
1913	1 2
1914	2 10
1915	3 104
1916	"""
1917	TE = []
1918	for i in range(n):
1919	for j in range(n):
1920	TE.append((i, j))
1921	m = len(TE)
1922	Glist= []
1923	for i in range(2**m):
1924	G = DiGraph(loops=True)
1925	G.add_vertices(range(n))
1926	b = Integer(i).binary()
1927	b = '0'*(m-len(b)) + b
1928	for j in range(m):
1929	if int(b[j]):
1930	G.add_edge(TE[j])
1931	Glist.append(G)
1932	return Glist
1933
1934	def all_labeled_digraphs(n):
1935	"""
1936	EXAMPLES:
1937	sage: import sage.graphs.graph_isom
1938	sage: from sage.graphs.graph_isom import search_tree, all_labeled_digraphs
1939	sage: from sage.graphs.graph import enum
1940	sage: Glist = {}
1941	sage: Giso = {}
1942	sage: for n in range(1,4):
1943	... Glist[n] = all_labeled_digraphs(n)
1944	... Giso[n] = []
1945	... for g in Glist[n]:
1946	... a, b = search_tree(g, [range(n)], dig=True)
1947	... inn = False
1948	... for gi in Giso[n]:
1949	... if enum(b) == enum(gi):
1950	... inn = True
1951	... if not inn:
1952	... Giso[n].append(b)
1953	sage: for n in Giso:
1954	... print n, len(Giso[n])
1955	1 1
1956	2 3
1957	3 16
1958	"""
1959
1960	## sage: n = 4 # long time (4 minutes)
1961	## sage: Glist[n] = all_labeled_digraphs(n) # long time
1962	## sage: Giso[n] = [] # long time
1963	## sage: for g in Glist[n]: # long time
1964	## ... a, b = search_tree(g, [range(n)], dig=True)
1965	## ... inn = False
1966	## ... for gi in Giso[n]:
1967	## ... if enum(b) == enum(gi):
1968	## ... inn = True
1969	## ... if not inn:
1970	## ... Giso[n].append(b)
1971	## sage: print n, len(Giso[n]) # long time
1972	## 4 218
1973	TE = []
1974	for i in range(n):
1975	for j in range(n):
1976	if i != j: TE.append((i, j))
1977	m = len(TE)
1978	Glist= []
1979	for i in range(2**m):
1980	G = DiGraph(loops=True)
1981	G.add_vertices(range(n))
1982	b = Integer(i).binary()
1983	b = '0'*(m-len(b)) + b
1984	for j in range(m):
1985	if int(b[j]):
1986	G.add_edge(TE[j])
1987	Glist.append(G)
1988	return Glist
1989
1990	def perm_group_elt(lperm):
1991	"""
1992	Given a list permutation of the set {0, 1, ..., n-1},
1993	returns the corresponding PermutationGroupElement where
1994	we take 0 = n.
1995	"""
1996	from sage.groups.perm_gps.permgroup_named import SymmetricGroup
1997	n = len(lperm)
1998	S = SymmetricGroup(n)
1999	Part = orbit_partition(lperm, list_perm=True)
2000	gens = []
2001	for z in Part:
2002	if len(z) > 1:
2003	if 0 in z:
2004	zed = z.index(0)
2005	generator = z[:zed] + [n] + z[zed+1:]
2006	gens.append(tuple(generator))
2007	else:
2008	gens.append(tuple(z))
2009	E = S(gens)
2010	return E
2011
2012	def orbit_partition(gamma, list_perm=False):
2013	r"""
2014	Assuming that G is a graph on vertices {0,1,...,n-1}, and gamma is an
2015	element of SymmetricGroup(n), returns the partition of the vertex set
2016	determined by the orbits of gamma, considered as action on the set
2017	{1,2,...,n} where we take 0 = n. In other words, returns the partition
2018	determined by a cyclic representation of gamma.
2019
2020	INPUT:
2021	list_perm -- if True, assumes gamma is a list representing the map
2022	i \mapsto gamma[i].
2023
2024	EXAMPLES:
2025	sage: import sage.graphs.graph_isom
2026	sage: from sage.graphs.graph_isom import orbit_partition
2027	sage: G = graphs.PetersenGraph()
2028	sage: S = SymmetricGroup(10)
2029	sage: gamma = S('(10,1,2,3,4)(5,6,7)(8,9)')
2030	sage: orbit_partition(gamma)
2031	[[1, 2, 3, 4, 0], [5, 6, 7], [8, 9]]
2032	sage: gamma = S('(10,5)(1,6)(2,7)(3,8)(4,9)')
2033	sage: orbit_partition(gamma)
2034	[[1, 6], [2, 7], [3, 8], [4, 9], [5, 0]]
2035	"""
2036	if list_perm:
2037	n = len(gamma)
2038	seen = [1] + [0]*(n-1)
2039	i = 0
2040	p = 0
2041	partition = [[0]]
2042	while sum(seen) < n:
2043	if gamma[i] != partition[p][0]:
2044	partition[p].append(gamma[i])
2045	i = gamma[i]
2046	seen[i] = 1
2047	else:
2048	i = min([j for j in range(n) if seen[j] == 0])
2049	partition.append([i])
2050	p += 1
2051	seen[i] = 1
2052	return partition
2053	else:
2054	n = len(gamma.list())
2055	l = []
2056	for i in range(1,n+1):
2057	orb = gamma.orbit(i)
2058	if orb not in l: l.append(orb)
2059	for i in l:
2060	for j in range(len(i)):
2061	if i[j] == n:
2062	i[j] = 0
2063	return l
2064
2065	def number_of_graphs(n, j = None):
2066	graph_list = []
2067	n = int(n)
2068	l = 2*((n(n-1))/2)
2069	print 'Computing canonical labels for %d labeled graphs.'%l
2070	k = 0
2071	l = l/10
2072	if l > 100: l = 100
2073	for g in all_labeled_graphs(n):
2074	if k%l == 0:
2075	print k
2076	k += 1
2077	g = g.canonical_label()
2078	if g not in graph_list:
2079	graph_list.append(g)
2080	return len(graph_list)
2081

Note: See TracBrowser for help on using the repository browser.

Download in other formats: