CoCalc -- refinement

GitHub Repository: sagemath/sagelib
Path: blob/master/sage/groups/perm_gps/partn_ref/refinement_binary.pyx
⁴⁰⁵⁷ views
1
"""
2
Partition backtrack functions for binary codes
3

4
DOCTEST:
5
    sage: import sage.groups.perm_gps.partn_ref.refinement_binary
6

7
REFERENCE:
8

9
    [1] McKay, Brendan D. Practical Graph Isomorphism. Congressus Numerantium,
10
        Vol. 30 (1981), pp. 45-87.
11

12
    [2] Leon, Jeffrey. Permutation Group Algorithms Based on Partitions, I:
13
        Theory and Algorithms. J. Symbolic Computation, Vol. 12 (1991), pp.
14
        533-583.
15

16
"""
17

18
#*****************************************************************************
19
#      Copyright (C) 2006 - 2011 Robert L. Miller <[email protected]>
20
#
21
# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
22
#                         http://www.gnu.org/licenses/
23
#*****************************************************************************
24

25
include 'data_structures_pyx.pxi' # includes bitsets
26

27
from sage.matrix.matrix import is_Matrix
28

29
cdef class LinearBinaryCodeStruct(BinaryCodeStruct):
30
    
31
    def __cinit__(self, matrix):
32
        cdef int i,j
33
        self.degree = matrix.ncols()
34
        self.dimension = matrix.nrows()
35
        if self.dimension >= (sizeof(int) << 3):
36
            raise NotImplementedError
37
            # By the time the dimension gets this big, the computation is infeasible anyway...
38
        self.nwords = 1<<self.dimension
39

40
        self.basis = <bitset_s *> sage_malloc(self.dimension * sizeof(bitset_s))
41
        self.scratch_bitsets = <bitset_s *> sage_malloc((2*self.dimension+2) * sizeof(bitset_s))
42
        self.alpha_is_wd = <bitset_s *> sage_malloc(sizeof(bitset_s))
43
        self.word_ps = PS_new(self.nwords, 1)
44
        self.alpha = <int *> sage_malloc((self.nwords+self.degree) * sizeof(int))
45
        self.scratch = <int *> sage_malloc((3*self.nwords+3*self.degree+2) * sizeof(int))
46
        
47
        if self.basis       is NULL or self.scratch_bitsets is NULL \
48
        or self.alpha_is_wd is NULL or self.word_ps         is NULL \
49
        or self.alpha       is NULL or self.scratch         is NULL:
50
            if self.basis is not NULL: sage_free(self.basis)
51
            if self.scratch_bitsets is not NULL: sage_free(self.scratch_bitsets)
52
            if self.alpha_is_wd is not NULL: sage_free(self.alpha_is_wd)
53
            if self.word_ps is not NULL: PS_dealloc(self.word_ps)
54
            if self.alpha is not NULL: sage_free(self.alpha)
55
            if self.scratch is not NULL: sage_free(self.scratch)
56
            raise MemoryError
57
        
58
        cdef bint memerr = 0
59
        for i from 0 <= i < self.dimension:
60
            try: bitset_init(&self.basis[i], self.degree)
61
            except MemoryError:
62
                for j from 0 <= j < i:
63
                    bitset_free(&self.basis[j])
64
                memerr = 1
65
                break
66
        if not memerr:
67
            for i from 0 <= i < 2*self.dimension+2:
68
                try: bitset_init(&self.scratch_bitsets[i], self.degree)
69
                except MemoryError:
70
                    for j from 0 <= j < i:
71
                        bitset_free(&self.scratch_bitsets[j])
72
                    for j from 0 <= j < self.dimension:
73
                        bitset_free(&self.basis[j])
74
                    memerr = 1
75
                    break
76
        if not memerr:
77
            try: bitset_init(self.alpha_is_wd, self.nwords + self.degree)
78
            except MemoryError:
79
                for j from 0 <= j < 2*self.dimension+2:
80
                    bitset_free(&self.scratch_bitsets[j])
81
                for j from 0 <= j < self.dimension:
82
                    bitset_free(&self.basis[j])
83
                memerr = 1
84
        if memerr:
85
            sage_free(self.basis); sage_free(self.scratch_bitsets)
86
            sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
87
            sage_free(self.alpha); sage_free(self.scratch)
88
            raise MemoryError
89
        else:
90
            bitset_zero(self.alpha_is_wd)
91
            for j from 0 <= j < self.dimension:
92
                bitset_zero(&self.basis[j])
93
        
94
        for i,j in matrix.nonzero_positions():
95
            bitset_set(&self.basis[i], j)
96
        
97
        self.output = NULL
98
        self.ith_word = &ith_word_linear
99
    
100
    def run(self, partition=None):
101
        """
102
        Perform the canonical labeling and automorphism group computation,
103
        storing results to self.
104
        
105
        INPUT:
106
        partition -- an optional list of lists partition of the columns.
107
            default is the unit partition.
108
        
109
        EXAMPLES:
110
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct
111
            
112
            sage: B = LinearBinaryCodeStruct(matrix(GF(2),[[1,0,1],[0,1,1]]))
113
            sage: B.run()
114
            sage: B.automorphism_group()
115
            ([[0, 2, 1], [1, 0, 2]], 6, [0, 1])
116
            sage: B.canonical_relabeling()
117
            [0, 1, 2]
118

119
            sage: B = LinearBinaryCodeStruct(matrix(GF(2),[[1,1,1,1]]))
120
            sage: B.automorphism_group()
121
            ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], 24, [0, 1, 2])
122
            sage: B.canonical_relabeling()
123
            [0, 1, 2, 3]
124

125
            sage: B = LinearBinaryCodeStruct(matrix(GF(2),[[0,0,0,0,0,0,0,0,0,0,0,0,0,0,1]]))
126
            sage: B.automorphism_group()[1] == factorial(14)
127
            True
128

129
            sage: M = Matrix(GF(2),[\
130
            ... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0],\
131
            ... [0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0],\
132
            ... [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1],\
133
            ... [0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],\
134
            ... [0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1]])
135
            sage: B = LinearBinaryCodeStruct(M)
136
            sage: B.automorphism_group()[1]
137
            322560
138

139
            sage: M = Matrix(GF(2),[\
140
            ... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0],\
141
            ... [0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0],\
142
            ... [0,0,0,0,0,1,0,1,0,0,0,1,1,1,1,1,1],\
143
            ... [0,0,0,1,1,0,0,0,0,1,1,0,1,1,0,1,1]])
144
            sage: B = LinearBinaryCodeStruct(M)
145
            sage: B.automorphism_group()[1]
146
            2304
147

148
            sage: M=Matrix(GF(2),[\
149
            ... [1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0,0],\
150
            ... [0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0,0],\
151
            ... [0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0,0],\
152
            ... [0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0,0],\
153
            ... [0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0,0],\
154
            ... [0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0,0],\
155
            ... [0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1,0],\
156
            ... [0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,1]])
157
            sage: B = LinearBinaryCodeStruct(M)
158
            sage: B.automorphism_group()[1]
159
            136
160

161
            sage: M = Matrix(GF(2),[\
162
            ... [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0],
163
            ... [0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0],
164
            ... [0,0,0,0,1,1,0,0,0,0,0,0,1,1,1,1,1,1],
165
            ... [0,0,1,1,0,0,0,0,0,0,1,1,1,1,0,0,1,1],
166
            ... [0,0,0,1,0,0,0,1,0,1,0,1,0,1,1,1,0,1],
167
            ... [0,1,0,0,0,1,0,0,0,1,1,1,0,1,0,1,1,0]])
168
            sage: B = LinearBinaryCodeStruct(M)
169
            sage: B.automorphism_group()[1]
170
            2160
171

172
            sage: M=Matrix(GF(2),[\
173
            ... [0,1,0,1,1,1,0,0,0,1,0,0,0,1,0,0,0,1,1,1,0,1],\
174
            ... [1,0,1,1,1,0,0,0,1,0,0,0,1,0,0,0,1,1,1,0,1,0],\
175
            ... [0,1,1,1,0,0,0,1,0,0,1,1,0,0,0,1,1,1,0,1,0,0],\
176
            ... [1,1,1,0,0,0,1,0,0,1,0,0,0,0,1,1,1,0,1,0,0,1],\
177
            ... [1,1,0,0,0,1,0,0,1,0,1,0,0,1,1,1,0,1,0,0,1,0],\
178
            ... [1,0,0,0,1,0,0,1,0,1,1,0,1,1,1,0,1,0,0,1,0,0],\
179
            ... [0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,0,1,0,0,0],\
180
            ... [0,0,1,0,0,1,0,1,1,1,0,1,1,0,1,0,0,1,0,0,0,1],\
181
            ... [0,1,0,0,1,0,1,1,1,0,0,1,0,1,0,0,1,0,0,0,1,1],\
182
            ... [1,0,0,1,0,1,1,1,0,0,0,0,1,0,0,1,0,0,0,1,1,1],\
183
            ... [0,0,1,0,1,1,1,0,0,0,1,1,0,0,1,0,0,0,1,1,1,0]])
184
            sage: B = LinearBinaryCodeStruct(M)
185
            sage: B.automorphism_group()[1]
186
            887040
187

188
            sage: M = Matrix(GF(2),[\
189
            ... [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
190
            ... [0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
191
            ... [0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0],
192
            ... [0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0],
193
            ... [0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0],
194
            ... [0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0],
195
            ... [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1],
196
            ... [1,0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,1,1,0,0],
197
            ... [1,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,0,0],
198
            ... [1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,0,1,0]])
199
            sage: B = LinearBinaryCodeStruct(M)
200
            sage: B.automorphism_group()[1]
201
            294912
202

203
            sage: M = Matrix(GF(2), [\
204
            ... [1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
205
            ... [0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
206
            ... [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0],
207
            ... [0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,0],
208
            ... [0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,1,1,1,1,0],
209
            ... [0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0],
210
            ... [0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1],
211
            ... [0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,1,0,0,1,1,1,0,1],
212
            ... [0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,1,1,1,1,0,0,0,1]])
213
            sage: B = LinearBinaryCodeStruct(M)
214
            sage: B.automorphism_group()[1]
215
            442368
216

217
            sage: M = Matrix(GF(2), [\
218
            ... [1,1,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0],\
219
            ... [1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,1]])
220
            sage: B = LinearBinaryCodeStruct(M)
221
            sage: B.automorphism_group()[1]
222
            17868913969917295853568000000
223
 
224
        """
225
        cdef int i, n = self.degree
226
        cdef PartitionStack *part
227
        if partition is None:
228
            part = PS_new(n, 1)
229
        else:
230
            part = PS_from_list(partition)
231
        if part is NULL:
232
            raise MemoryError
233

234
        self.first_time = 1
235
        
236
        self.output = get_aut_gp_and_can_lab(<void *>self, part, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, 1, NULL)
237
        
238
        PS_dealloc(part)
239
    
240
    def automorphism_group(self):
241
        """
242
        Returns a list of generators of the automorphism group, along with its
243
        order and a base for which the list of generators is a strong generating
244
        set.
245
        
246
        EXAMPLE: (For more examples, see self.run())
247
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct
248

249
            sage: B = LinearBinaryCodeStruct(matrix(GF(2),[[1,1,1,1]]))
250
            sage: B.automorphism_group()
251
            ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], 24, [0, 1, 2])
252

253
        """
254
        cdef int i, j
255
        cdef list generators, base
256
        cdef Integer order
257
        if self.output is NULL:
258
            self.run()
259
        generators = []
260
        for i from 0 <= i < self.output.num_gens:
261
            generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree])
262
        order = Integer()
263
        SC_order(self.output.group, 0, order.value)
264
        base = [self.output.group.base_orbits[i][0] for i from 0 <= i < self.output.group.base_size]
265
        return generators, order, base
266
    
267
    def canonical_relabeling(self):
268
        """
269
        Returns a canonical relabeling (in list permutation format).
270
        
271
        EXAMPLES: (For more examples, see self.run())
272
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct
273

274
            sage: B = LinearBinaryCodeStruct(matrix(GF(2), [[1,1,0]]))
275
            sage: B.automorphism_group()
276
            ([[1, 0, 2]], 2, [0])
277
            sage: B.canonical_relabeling()
278
            [1, 2, 0]
279
            sage: B = LinearBinaryCodeStruct(matrix(GF(2), [[1,0,1]]))
280
            sage: B.automorphism_group()
281
            ([[2, 1, 0]], 2, [0])
282
            sage: B.canonical_relabeling()
283
            [1, 0, 2]
284

285
        """
286
        cdef int i
287
        if self.output is NULL:
288
            self.run()
289
        return [self.output.relabeling[i] for i from 0 <= i < self.degree]
290
    
291
    def is_isomorphic(self, LinearBinaryCodeStruct other):
292
        """
293
        Calculate whether self is isomorphic to other.
294
        
295
        EXAMPLES:
296
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import LinearBinaryCodeStruct
297

298
            sage: B = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
299
            sage: C = LinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,1],[1,1,0,1,1,0]]))
300
            sage: B.is_isomorphic(C)
301
            [0, 1, 2, 5, 3, 4]
302
            
303
        """
304
        cdef int i, n = self.degree
305
        cdef int *output, *ordering
306
        cdef PartitionStack *part
307
        part = PS_new(n, 1)
308
        ordering = <int *> sage_malloc(self.degree * sizeof(int))
309
        if part is NULL or ordering is NULL:
310
            if part is not NULL: PS_dealloc(part)
311
            if ordering is not NULL: sage_free(ordering)
312
            raise MemoryError
313
        for i from 0 <= i < n:
314
            ordering[i] = i
315
        self.first_time = 1
316
        other.first_time = 1
317
        
318
        output = double_coset(<void *> self, <void *> other, part, ordering, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, NULL)
319
        
320
        PS_dealloc(part)
321
        sage_free(ordering)
322

323
        if output is NULL:
324
            return False
325
        else:
326
            output_py = [output[i] for i from 0 <= i < n]
327
            sage_free(output)
328
            return output_py
329
    
330
    def __dealloc__(self):
331
        cdef int j
332
        bitset_free(self.alpha_is_wd)
333
        for j from 0 <= j < 2*self.dimension+2:
334
            bitset_free(&self.scratch_bitsets[j])
335
        for j from 0 <= j < self.dimension:
336
            bitset_free(&self.basis[j])
337
        sage_free(self.basis); sage_free(self.scratch_bitsets)
338
        sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
339
        sage_free(self.alpha); sage_free(self.scratch)
340
        if self.output is not NULL:
341
            sage_free(self.output.generators)
342
            SC_dealloc(self.output.group)
343
            sage_free(self.output.relabeling)
344
            sage_free(self.output)
345
            
346
cdef int ith_word_linear(BinaryCodeStruct self, int i, bitset_s *word):
347
    cdef LinearBinaryCodeStruct LBCS = <LinearBinaryCodeStruct> self
348
    cdef int j
349
    bitset_zero(word)
350
    for j from 0 <= j < LBCS.dimension:
351
        if (1<<j)&i:
352
            bitset_xor(word, word, &LBCS.basis[j])
353
    return 0
354

355
cdef class NonlinearBinaryCodeStruct(BinaryCodeStruct):
356
    
357
    def __cinit__(self, arg):
358
        cdef int i,j
359
        if is_Matrix(arg):
360
            self.degree = arg.ncols()
361
            self.nwords = arg.nrows()
362
        elif isinstance(arg, tuple):
363
            assert len(arg) == 2
364
            self.degree, self.nwords = arg
365
        else:
366
            raise NotImplementedError
367

368
        self.words = <bitset_s *> sage_malloc(self.nwords * sizeof(bitset_s))
369
        self.scratch_bitsets = <bitset_s *> sage_malloc((4*self.nwords+1) * sizeof(bitset_s))
370
        self.alpha_is_wd = <bitset_s *> sage_malloc(sizeof(bitset_s))
371
        self.word_ps = PS_new(self.nwords, 1)
372
        self.alpha = <int *> sage_malloc((self.nwords+self.degree) * sizeof(int))
373
        self.scratch = <int *> sage_malloc((3*self.nwords+3*self.degree+2) * sizeof(int))
374
        if self.words       is NULL or self.scratch_bitsets is NULL \
375
        or self.alpha_is_wd is NULL or self.word_ps         is NULL \
376
        or self.alpha       is NULL or self.scratch         is NULL:
377
            if self.words is not NULL: sage_free(self.words)
378
            if self.scratch_bitsets is not NULL: sage_free(self.scratch_bitsets)
379
            if self.alpha_is_wd is not NULL: sage_free(self.alpha_is_wd)
380
            if self.word_ps is not NULL: PS_dealloc(self.word_ps)
381
            if self.alpha is not NULL: sage_free(self.alpha)
382
            if self.scratch is not NULL: sage_free(self.scratch)
383
            raise MemoryError
384
        
385
        cdef bint memerr = 0
386
        for i from 0 <= i < self.nwords:
387
            try: bitset_init(&self.words[i], self.degree)
388
            except MemoryError:
389
                for j from 0 <= j < i:
390
                    bitset_free(&self.words[j])
391
                memerr = 1
392
                break
393
        if not memerr:
394
            for i from 0 <= i < 4*self.nwords:
395
                try: bitset_init(&self.scratch_bitsets[i], self.degree)
396
                except MemoryError:
397
                    for j from 0 <= j < i:
398
                        bitset_free(&self.scratch_bitsets[j])
399
                    for j from 0 <= j < self.nwords:
400
                        bitset_free(&self.words[j])
401
                    memerr = 1
402
                    break
403
        if not memerr:
404
            try: bitset_init(&self.scratch_bitsets[4*self.nwords], self.nwords)
405
            except MemoryError:
406
                for j from 0 <= j < 4*self.nwords:
407
                    bitset_free(&self.scratch_bitsets[j])
408
                for j from 0 <= j < self.nwords:
409
                    bitset_free(&self.words[j])
410
                memerr = 1
411
        if not memerr:
412
            try: bitset_init(self.alpha_is_wd, self.nwords + self.degree)
413
            except MemoryError:
414
                for j from 0 <= j < 4*self.nwords + 1:
415
                    bitset_free(&self.scratch_bitsets[j])
416
                for j from 0 <= j < self.nwords:
417
                    bitset_free(&self.words[j])
418
                memerr = 1
419
        if memerr:
420
            sage_free(self.words); sage_free(self.scratch_bitsets)
421
            sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
422
            sage_free(self.alpha); sage_free(self.scratch)
423
            raise MemoryError
424
        else:
425
            bitset_zero(self.alpha_is_wd)
426
            for j from 0 <= j < self.nwords:
427
                bitset_zero(&self.words[j])
428
        
429
        if is_Matrix(arg):
430
            for i,j in arg.nonzero_positions():
431
                bitset_set(&self.words[i], j)
432
        
433
        self.output = NULL
434
        self.ith_word = &ith_word_nonlinear
435

436
    def __dealloc__(self):
437
        cdef int j
438
        bitset_free(self.alpha_is_wd)
439
        for j from 0 <= j < 4*self.nwords + 1:
440
            bitset_free(&self.scratch_bitsets[j])
441
        for j from 0 <= j < self.nwords:
442
            bitset_free(&self.words[j])
443
        sage_free(self.words); sage_free(self.scratch_bitsets)
444
        sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
445
        sage_free(self.alpha); sage_free(self.scratch)
446
        if self.output is not NULL:
447
            sage_free(self.output.generators)
448
            SC_dealloc(self.output.group)
449
            sage_free(self.output.relabeling)
450
            sage_free(self.output)
451
            
452
    def run(self, partition=None):
453
        """
454
        Perform the canonical labeling and automorphism group computation,
455
        storing results to self.
456
        
457
        INPUT:
458
        partition -- an optional list of lists partition of the columns.
459
            default is the unit partition.
460
        
461
        EXAMPLES:
462
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
463

464
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,0,0,0],[0,0,1,0]]))
465
            sage: B.run()
466
            sage: B.automorphism_group()
467
            ([[2, 1, 0, 3], [0, 3, 2, 1]], 4, [1, 0])
468
            sage: B.canonical_relabeling()
469
            [2, 0, 3, 1]
470

471
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0],[1,1,0,1],[1,0,1,1],[0,1,1,1]]))
472
            sage: B.run()
473
            sage: B.automorphism_group()
474
            ([[0, 1, 3, 2], [0, 2, 1, 3], [1, 0, 2, 3]], 24, [0, 1, 2])
475
            sage: B.canonical_relabeling()
476
            [0, 1, 2, 3]
477

478
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,0],[1,1,0,1,0,0],[1,0,1,1,0,0],[0,1,1,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1]]))
479
            sage: B.run()
480
            sage: B.automorphism_group()
481
            ([[0, 1, 3, 2, 4, 5],
482
              [0, 2, 1, 3, 4, 5],
483
              [1, 0, 2, 3, 4, 5],
484
              [0, 1, 2, 3, 5, 4]],
485
             48,
486
             [4, 0, 1, 2])
487
            sage: B.canonical_relabeling()
488
            [2, 3, 4, 5, 0, 1]
489

490
        """
491
        cdef int n = self.degree
492
        cdef PartitionStack *part
493
        if partition is None:
494
            part = PS_new(n, 1)
495
        else:
496
            part = PS_from_list(partition)
497
        if part is NULL:
498
            raise MemoryError
499
        self.first_time = 1
500
        
501
        self.output = get_aut_gp_and_can_lab(<void *> self, part, self.degree, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, 1, NULL)
502
        
503
        PS_dealloc(part)
504
    
505
    def automorphism_group(self):
506
        """
507
        Returns a list of generators of the automorphism group, along with its
508
        order and a base for which the list of generators is a strong generating
509
        set.
510
        
511
        EXAMPLE: (For more examples, see self.run())
512
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
513

514
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,0],[1,1,0,1,0,0],[1,0,1,1,0,0],[0,1,1,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1]]))
515
            sage: B.run()
516
            sage: B.automorphism_group()
517
            ([[0, 1, 3, 2, 4, 5],
518
              [0, 2, 1, 3, 4, 5],
519
              [1, 0, 2, 3, 4, 5],
520
              [0, 1, 2, 3, 5, 4]],
521
             48,
522
             [4, 0, 1, 2])
523

524
        """
525
        cdef int i, j
526
        cdef list generators, base
527
        cdef Integer order
528
        if self.output is NULL:
529
            self.run()
530
        generators = []
531
        for i from 0 <= i < self.output.num_gens:
532
            generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree])
533
        order = Integer()
534
        SC_order(self.output.group, 0, order.value)
535
        base = [self.output.group.base_orbits[i][0] for i from 0 <= i < self.output.group.base_size]
536
        return generators, order, base
537
    
538
    def canonical_relabeling(self):
539
        """
540
        Returns a canonical relabeling (in list permutation format).
541
        
542
        EXAMPLES: (For more examples, see self.run())
543
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
544

545
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,0,0,0],[1,1,0,1,0,0],[1,0,1,1,0,0],[0,1,1,1,0,0],[0,0,0,0,1,0],[0,0,0,0,0,1]]))
546
            sage: B.run()
547
            sage: B.canonical_relabeling()
548
            [2, 3, 4, 5, 0, 1]
549

550
        """
551
        cdef int i
552
        if self.output is NULL:
553
            self.run()
554
        return [self.output.relabeling[i] for i from 0 <= i < self.degree]
555

556
    def is_isomorphic(self, NonlinearBinaryCodeStruct other):
557
        """
558
        Calculate whether self is isomorphic to other.
559
        
560
        EXAMPLES:
561
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
562

563
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
564
            sage: C = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,0,0,1,1],[1,1,1,1,0,0]]))
565
            sage: B.is_isomorphic(C)
566
            [2, 3, 0, 1, 4, 5]
567
            
568
        """
569
        cdef int i, n = self.degree
570
        cdef int *output, *ordering
571
        cdef PartitionStack *part
572
        part = PS_new(n, 1)
573
        ordering = <int *> sage_malloc(n * sizeof(int))
574
        if part is NULL or ordering is NULL:
575
            if part is not NULL: PS_dealloc(part)
576
            if ordering is not NULL: sage_free(ordering)
577
            raise MemoryError
578
        for i from 0 <= i < n:
579
            ordering[i] = i
580
        self.first_time = 1
581
        other.first_time = 1
582
        
583
        output = double_coset(<void *> self, <void *> other, part, ordering, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, NULL)
584
        
585
        PS_dealloc(part)
586
        sage_free(ordering)
587

588
        if output is NULL:
589
            return False
590
        else:
591
            output_py = [output[i] for i from 0 <= i < self.degree]
592
            sage_free(output)
593
            return output_py
594
    
595
cdef int ith_word_nonlinear(BinaryCodeStruct self, int i, bitset_s *word):
596
    cdef NonlinearBinaryCodeStruct NBCS = <NonlinearBinaryCodeStruct> self
597
    bitset_copy(word, &NBCS.words[i])
598
    return 0
599

600
cdef int refine_by_bip_degree(PartitionStack *col_ps, void *S, int *cells_to_refine_by, int ctrb_len):
601
    """
602
    Refines the input partition by checking degrees of vertices to the given
603
    cells in the associated bipartite graph (vertices split into columns and
604
    words).
605
    
606
    INPUT:
607
    col_ps -- a partition stack, whose finest partition is the partition to be
608
        refined.
609
    S -- a binary code struct object
610
    cells_to_refine_by -- a list of pointers to cells to check degrees against
611
        in refining the other cells (updated in place)
612
    ctrb_len -- how many cells in cells_to_refine_by
613
    
614
    OUTPUT:
615

616
    An integer $I$ invariant under the orbits of $S_n$.  That is, if $\gamma$ is a 
617
    permutation of the columns, then 
618
    $$ I(G, PS, cells_to_refine_by) = I( \gamma(G), \gamma(PS), \gamma(cells_to_refine_by) ) .$$
619
    
620
    """
621
    cdef BinaryCodeStruct BCS = <BinaryCodeStruct> S
622
    cdef int current_cell_against = 0
623
    cdef int current_cell, i, r, j
624
    cdef int first_largest_subcell
625
    cdef int invariant = 0
626
    cdef PartitionStack *word_ps = BCS.word_ps
627
    cdef int *ctrb = BCS.alpha
628
    cdef bitset_s *ctrb_is_wd = BCS.alpha_is_wd
629

630
    word_ps.depth = col_ps.depth
631
    PS_clear(word_ps)
632
    bitset_zero(ctrb_is_wd)
633
    memcpy(ctrb, cells_to_refine_by, ctrb_len * sizeof(int))
634
    if BCS.first_time:
635
        BCS.first_time = 0
636
        ctrb[ctrb_len] = 0
637
        bitset_set(ctrb_is_wd, ctrb_len)
638
        ctrb_len += 1
639
    cdef int *col_degrees = BCS.scratch                                     # len degree
640
    cdef int *col_counts  = &BCS.scratch[BCS.degree]                        # len nwords+1
641
    cdef int *col_output  = &BCS.scratch[BCS.degree + BCS.nwords + 1]       # len degree
642
    cdef int *word_degrees  = &BCS.scratch[2*BCS.degree + BCS.nwords + 1]   # len nwords
643
    cdef int *word_counts   = &BCS.scratch[2*BCS.degree + 2*BCS.nwords + 1] # len degree+1
644
    cdef int *word_output   = &BCS.scratch[3*BCS.degree + 2*BCS.nwords + 2] # len nwords
645
    cdef bint necessary_to_split_cell
646
    cdef int against_index
647
    while not (PS_is_discrete(col_ps) and PS_is_discrete(word_ps)) and current_cell_against < ctrb_len:
648
        invariant += 1
649
        current_cell = 0
650
        if bitset_check(ctrb_is_wd, current_cell_against):
651
            while current_cell < col_ps.degree:
652
                invariant += 8
653
                i = current_cell
654
                necessary_to_split_cell = 0
655
                while 1:
656
                    col_degrees[i-current_cell] = col_degree(col_ps, BCS, i, ctrb[current_cell_against], word_ps)
657
                    if col_degrees[i-current_cell] != col_degrees[0]:
658
                        necessary_to_split_cell = 1
659
                    i += 1
660
                    if col_ps.levels[i-1] <= col_ps.depth:
661
                        break
662
                # now, i points to the next cell (before refinement)
663
                if necessary_to_split_cell:
664
                    invariant += 8
665
                    first_largest_subcell = sort_by_function_codes(col_ps, current_cell, col_degrees, col_counts, col_output, BCS.nwords+1)
666
                    invariant += col_degree(col_ps, BCS, i-1, ctrb[current_cell_against], word_ps)
667
                    invariant += first_largest_subcell
668
                    against_index = current_cell_against
669
                    while against_index < ctrb_len:
670
                        if ctrb[against_index] == current_cell and not bitset_check(ctrb_is_wd, against_index):
671
                            ctrb[against_index] = first_largest_subcell
672
                            break
673
                        against_index += 1
674
                    r = current_cell
675
                    while 1:
676
                        if r == current_cell or col_ps.levels[r-1] == col_ps.depth:
677
                            if r != first_largest_subcell:
678
                                ctrb[ctrb_len] = r
679
                                ctrb_len += 1
680
                        r += 1
681
                        if r >= i:
682
                            break
683
                    invariant += (i - current_cell)
684
                current_cell = i
685
        else:
686
            while current_cell < word_ps.degree:
687
                invariant += 64
688
                i = current_cell
689
                necessary_to_split_cell = 0
690
                while 1:
691
                    word_degrees[i-current_cell] = word_degree(word_ps, BCS, i, ctrb[current_cell_against], col_ps)
692
                    if word_degrees[i-current_cell] != word_degrees[0]:
693
                        necessary_to_split_cell = 1
694
                    i += 1
695
                    if word_ps.levels[i-1] <= col_ps.depth:
696
                        break
697
                # now, i points to the next cell (before refinement)
698
                if necessary_to_split_cell:
699
                    invariant += 64
700
                    first_largest_subcell = sort_by_function_codes(word_ps, current_cell, word_degrees, word_counts, word_output, BCS.degree+1)
701
                    invariant += word_degree(word_ps, BCS, i-1, ctrb[current_cell_against], col_ps)
702
                    invariant += first_largest_subcell
703
                    against_index = current_cell_against
704
                    while against_index < ctrb_len:
705
                        if ctrb[against_index] == current_cell and bitset_check(ctrb_is_wd, against_index):
706
                            ctrb[against_index] = first_largest_subcell
707
                            break
708
                        against_index += 1
709
                    r = current_cell
710
                    while 1:
711
                        if r == current_cell or word_ps.levels[r-1] == col_ps.depth:
712
                            if r != first_largest_subcell:
713
                                ctrb[ctrb_len] = r
714
                                bitset_set(ctrb_is_wd, ctrb_len)
715
                                ctrb_len += 1
716
                        r += 1
717
                        if r >= i:
718
                            break
719
                    invariant += (i - current_cell)
720
                current_cell = i
721
        current_cell_against += 1
722
    return invariant
723

724
cdef int compare_linear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2):
725
    """
726
    Compare gamma_1(S1) and gamma_2(S2).
727
    
728
    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
729
    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
730
    
731
    Abstractly, what this function does is relabel the basis of B by gamma_1 and
732
    gamma_2, run a row reduction on each, and verify that the matrices are the
733
    same, which holds if and only if the rowspan is the same. In practice, if
734
    the codes are not equal, the reductions (which are carried out in an
735
    interleaved way) will terminate as soon as this is discovered, and whichever
736
    code has a 1 in the entry in which they differ is reported as larger.
737
    
738
    INPUT:
739
    gamma_1, gamma_2 -- list permutations (inverse)
740
    S1, S2 -- binary code struct objects
741
    
742
    """
743
    cdef int i, piv_loc_1, piv_loc_2, cur_col, cur_row=0
744
    cdef bint is_pivot_1, is_pivot_2
745
    cdef LinearBinaryCodeStruct BCS1 = <LinearBinaryCodeStruct> S1
746
    cdef LinearBinaryCodeStruct BCS2 = <LinearBinaryCodeStruct> S2
747
    cdef bitset_s *basis_1 = BCS1.scratch_bitsets                   # len = dim
748
    cdef bitset_s *basis_2 = &BCS1.scratch_bitsets[BCS1.dimension]   # len = dim
749
    cdef bitset_s *pivots  = &BCS1.scratch_bitsets[2*BCS1.dimension] # len 1
750
    cdef bitset_s *temp  = &BCS1.scratch_bitsets[2*BCS1.dimension+1] # len 1
751
    for i from 0 <= i < BCS1.dimension:
752
        bitset_copy(&basis_1[i], &BCS1.basis[i])
753
        bitset_copy(&basis_2[i], &BCS2.basis[i])
754
    bitset_zero(pivots)
755
    for cur_col from 0 <= cur_col < BCS1.degree:
756
        is_pivot_1 = 0
757
        is_pivot_2 = 0
758
        for i from cur_row <= i < BCS1.dimension:
759
            if bitset_check(&basis_1[i], gamma_1[cur_col]):
760
                is_pivot_1 = 1
761
                piv_loc_1 = i
762
                break
763
        for i from cur_row <= i < BCS1.dimension:
764
            if bitset_check(&basis_2[i], gamma_2[cur_col]):
765
                is_pivot_2 = 1
766
                piv_loc_2 = i
767
                break
768
        if is_pivot_1 != is_pivot_2:
769
            return <int>is_pivot_2 - <int>is_pivot_1
770
        if is_pivot_1:
771
            bitset_set(pivots, cur_col)
772
            if piv_loc_1 != cur_row:
773
                bitset_copy(temp, &basis_1[piv_loc_1])
774
                bitset_copy(&basis_1[piv_loc_1], &basis_1[cur_row])
775
                bitset_copy(&basis_1[cur_row], temp)
776
            if piv_loc_2 != cur_row:
777
                bitset_copy(temp, &basis_2[piv_loc_2])
778
                bitset_copy(&basis_2[piv_loc_2], &basis_2[cur_row])
779
                bitset_copy(&basis_2[cur_row], temp)
780
            for i from 0 <= i < cur_row:
781
                if bitset_check(&basis_1[i], gamma_1[cur_col]):
782
                    bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
783
                if bitset_check(&basis_2[i], gamma_2[cur_col]):
784
                    bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row])
785
            for i from cur_row < i < BCS1.dimension:
786
                if bitset_check(&basis_1[i], gamma_1[cur_col]):
787
                    bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
788
                if bitset_check(&basis_2[i], gamma_2[cur_col]):
789
                    bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row])
790
            cur_row += 1
791
        else:
792
            for i from 0 <= i < cur_row:
793
                if bitset_check(&basis_1[i], gamma_1[cur_col]) != bitset_check(&basis_2[i], gamma_2[cur_col]):
794
                    return <int>bitset_check(&basis_2[i], gamma_2[cur_col]) - <int>bitset_check(&basis_1[i], gamma_1[cur_col])
795
    return 0
796

797
cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2):
798
    """
799
    Compare gamma_1(S1) and gamma_2(S2).
800
    
801
    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
802
    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
803
    
804
    INPUT:
805
    gamma_1, gamma_2 -- list permutations (inverse)
806
    S1, S2 -- a binary code struct object
807
    
808
    """
809
    cdef int side=0, i, start, end, n_one_1, n_one_2, cur_col
810
    cdef int where_0, where_1
811
    cdef NonlinearBinaryCodeStruct BCS1 = <NonlinearBinaryCodeStruct> S1
812
    cdef NonlinearBinaryCodeStruct BCS2 = <NonlinearBinaryCodeStruct> S2
813
    cdef bitset_s *B_1_0 = BCS1.scratch_bitsets                   # nwords of len degree
814
    cdef bitset_s *B_1_1 = &BCS1.scratch_bitsets[BCS1.nwords]      # nwords of len degree
815
    cdef bitset_s *B_2_0 = &BCS1.scratch_bitsets[2*BCS1.nwords]    # nwords of len degree
816
    cdef bitset_s *B_2_1 = &BCS1.scratch_bitsets[3*BCS1.nwords]    # nwords of len degree
817
    cdef bitset_s *dividers = &BCS1.scratch_bitsets[4*BCS1.nwords] # 1 of len nwords
818
    cdef bitset_s *B_1_this, *B_1_other, *B_2_this, *B_2_other
819
    for i from 0 <= i < BCS1.nwords:
820
        bitset_copy(&B_1_0[i], &BCS1.words[i])
821
        bitset_copy(&B_2_0[i], &BCS2.words[i])
822
    bitset_zero(dividers)
823
    bitset_set(dividers, BCS1.nwords-1)
824
    
825
    for cur_col from 0 <= cur_col < BCS1.degree:
826
        if side == 0:
827
            B_1_this  = B_1_0
828
            B_1_other = B_1_1
829
            B_2_this  = B_2_0
830
            B_2_other = B_2_1
831
        else:
832
            B_1_this  = B_1_1
833
            B_1_other = B_1_0
834
            B_2_this  = B_2_1
835
            B_2_other = B_2_0
836
        side ^= 1
837
        start = 0
838
        while start < BCS1.nwords:
839
            end = start
840
            while not bitset_check(dividers, end):
841
                end += 1
842
            end += 1
843
            n_one_1 = 0
844
            n_one_2 = 0
845
            for i from start <= i < end:
846
                n_one_1 += bitset_check(&B_1_this[i], gamma_1[cur_col])
847
                n_one_2 += bitset_check(&B_2_this[i], gamma_2[cur_col])
848
            if n_one_1 != n_one_2:
849
                if n_one_1 > n_one_2:
850
                    return 1
851
                else:
852
                    return -1
853
            where_0 = start
854
            where_1 = end - n_one_1
855
            if start < where_1 and where_1 < end:
856
                bitset_set(dividers, where_1 - 1)
857
            for i from start <= i < end:
858
                if bitset_check(&B_1_this[i], gamma_1[cur_col]):
859
                    bitset_copy(&B_1_other[where_1], &B_1_this[i])
860
                    where_1 += 1
861
                else:
862
                    bitset_copy(&B_1_other[where_0], &B_1_this[i])
863
                    where_0 += 1
864
            where_0 = start
865
            where_1 = end - n_one_2
866
            for i from start <= i < end:
867
                if bitset_check(&B_2_this[i], gamma_2[cur_col]):
868
                    bitset_copy(&B_2_other[where_1], &B_2_this[i])
869
                    where_1 += 1
870
                else:
871
                    bitset_copy(&B_2_other[where_0], &B_2_this[i])
872
                    where_0 += 1
873
            start = end
874

875
    return 0
876

877
cdef bint all_children_are_equivalent(PartitionStack *col_ps, void *S):
878
    """
879
    Returns True if any refinement of the current partition results in the same
880
    structure.
881
    
882
    WARNING:
883
    Converse does not hold in general! See Lemma 2.25 of [1] for details, noting
884
    that the binary code is interpreted as a bipartite graph (see module docs
885
    for details).
886
    
887
    INPUT:
888
    col_ps -- the partition stack to be checked
889
    S -- a binary code struct object
890
    """
891
    cdef BinaryCodeStruct BCS = <BinaryCodeStruct> S
892
    cdef PartitionStack *word_ps = BCS.word_ps
893
    cdef int i, n = col_ps.degree + BCS.nwords
894
    cdef bint in_cell = 0
895
    cdef int nontrivial_cells = 0
896
    cdef int total_cells = PS_num_cells(col_ps) + PS_num_cells(word_ps)
897
    if n <= total_cells + 4:
898
        return 1
899
    for i from 0 <= i < BCS.nwords:
900
        if word_ps.levels[i] <= col_ps.depth:
901
            if in_cell:
902
                nontrivial_cells += 1
903
            in_cell = 0
904
        else:
905
            in_cell = 1
906
    in_cell = 0
907
    for i from 0 <= i < BCS.degree:
908
        if col_ps.levels[i] <= col_ps.depth:
909
            if in_cell:
910
                nontrivial_cells += 1
911
            in_cell = 0
912
        else:
913
            in_cell = 1
914
    if n == total_cells + nontrivial_cells:
915
        return 1
916
    if n == total_cells + nontrivial_cells + 1:
917
        return 1
918
    return 0
919

920
cdef inline int word_degree(PartitionStack *word_ps, BinaryCodeStruct BCS, int entry, int cell_index, PartitionStack *col_ps):
921
    """
922
    Returns the number of edges from the vertex corresponding to entry to
923
    vertices in the cell corresponding to cell_index.
924
    
925
    INPUT:
926
    word_ps -- the partition stack to be checked
927
    col_ps -- corresponding partition stack on columns
928
    BCS -- a binary code struct object
929
    entry -- the position of the vertex in question in the entries of word_ps
930
    cell_index -- the starting position of the cell in question in the entries
931
        of PS
932
    """
933
    cdef bitset_t cell, word
934
    cdef int i, h
935
    bitset_init(cell, BCS.degree)
936
    bitset_zero(cell)
937
    bitset_init(word, BCS.degree)
938
    entry = word_ps.entries[entry]
939
    bitset_set(cell, col_ps.entries[cell_index])
940
    while col_ps.levels[cell_index] > col_ps.depth:
941
        cell_index += 1
942
        bitset_set(cell, col_ps.entries[cell_index])
943
    BCS.ith_word(BCS, entry, word)
944
    bitset_and(cell, word, cell)
945
    h = bitset_hamming_weight(cell)
946
    bitset_free(cell)
947
    bitset_free(word)
948
    return h
949

950
cdef inline int col_degree(PartitionStack *col_ps, BinaryCodeStruct BCS, int entry, int cell_index, PartitionStack *word_ps):
951
    """
952
    Returns the number of edges from the vertex corresponding to entry to
953
    vertices in the cell corresponding to cell_index.
954
    
955
    INPUT:
956
    col_ps -- the partition stack to be checked
957
    word_ps -- corresponding partition stack on words
958
    BCS -- a binary code struct object
959
    entry -- the position of the vertex in question in the entries of word_ps
960
    cell_index -- the starting position of the cell in question in the entries
961
        of PS
962
    """
963
    cdef bitset_t word
964
    bitset_init(word, BCS.degree)
965
    cdef int degree = 0, word_basis, i, b
966
    entry = col_ps.entries[entry]
967
    while 1:
968
        BCS.ith_word(BCS, word_ps.entries[cell_index], word)
969
        degree += bitset_check(word, entry)
970
        if not word_ps.levels[cell_index] > col_ps.depth:
971
            break
972
        cell_index += 1
973
    bitset_free(word)
974
    return degree
975

976
cdef inline int sort_by_function_codes(PartitionStack *PS, int start, int *degrees, int *counts, int *output, int count_max):
977
    """
978
    A simple counting sort, given the degrees of vertices to a certain cell.
979
    
980
    INPUT:
981
    PS -- the partition stack to be checked
982
    start -- beginning index of the cell to be sorted
983
    degrees -- the values to be sorted by
984
    count, count_max, output -- scratch space
985

986
    """
987
    cdef int n = PS.degree
988
    cdef int i, j, max, max_location
989
    for j from 0 <= j < count_max:
990
        counts[j] = 0
991
    i = 0
992
    while PS.levels[i+start] > PS.depth:
993
        counts[degrees[i]] += 1
994
        i += 1
995
    counts[degrees[i]] += 1
996
    # i+start is the right endpoint of the cell now
997
    max = counts[0]
998
    max_location = 0
999
    for j from 0 < j < count_max:
1000
        if counts[j] > max:
1001
            max = counts[j]
1002
            max_location = j
1003
        counts[j] += counts[j - 1]
1004
    for j from i >= j >= 0:
1005
        counts[degrees[j]] -= 1
1006
        output[counts[degrees[j]]] = PS.entries[start+j]
1007
    max_location = counts[max_location]+start
1008
    for j from 0 <= j <= i:
1009
        PS.entries[start+j] = output[j]
1010
    j = 1
1011
    while j < count_max and counts[j] <= i:
1012
        if counts[j] > 0:
1013
            PS.levels[start + counts[j] - 1] = PS.depth
1014
        PS_move_min_to_front(PS, start + counts[j-1], start + counts[j] - 1)
1015
        j += 1
1016
    return max_location
1017

1018

1019
def random_tests(num=50, n_max=50, k_max=6, nwords_max=200, perms_per_code=10, density_range=(.1,.9)):
1020
    """
1021
    Tests to make sure that C(gamma(B)) == C(B) for random permutations gamma
1022
    and random codes B, and that is_isomorphic returns an isomorphism.
1023

1024
    INPUT: 
1025
    num -- run tests for this many codes
1026
    n_max -- test codes with at most this many columns
1027
    k_max -- test codes with at most this for dimension
1028
    perms_per_code -- test each code with this many random permutations 
1029

1030
    DISCUSSION:
1031

1032
    This code generates num random linear codes B on at most n_max columns with
1033
    dimension at most k_max, and a random nonlinear code B2 on at most n_max
1034
    columns with number of words at most nwords_max. The density of entries in
1035
    the basis is chosen randomly between 0 and 1.
1036

1037
    For each code B (B2) generated, we uniformly generate perms_per_code random 
1038
    permutations and verify that the canonical labels of B and the image of B
1039
    under the generated permutation are equal, and check that the double coset
1040
    function returns an isomorphism.
1041
    
1042
    DOCTEST:
1043
        sage: import sage.groups.perm_gps.partn_ref.refinement_binary
1044
        sage: sage.groups.perm_gps.partn_ref.refinement_binary.random_tests()
1045
        All passed: ... random tests on ... codes.
1046

1047
    """
1048
    from sage.misc.misc import walltime
1049
    from sage.misc.prandom import random, randint
1050
    from sage.combinat.permutation import Permutations
1051
    from sage.matrix.constructor import random_matrix, matrix
1052
    from sage.rings.finite_rings.constructor import FiniteField as GF
1053
    cdef int h, i, j, n, k, num_tests = 0, num_codes = 0
1054
    cdef LinearBinaryCodeStruct B, C
1055
    cdef NonlinearBinaryCodeStruct B_n, C_n
1056
    for mmm in range(num):
1057
        p = random()*(density_range[1]-density_range[0]) + density_range[0]
1058
        n = randint(2, n_max)
1059
        k = randint(1, min(n-1,k_max) )
1060
        nwords = randint(1, min(n-1,nwords_max) )
1061
        S = Permutations(n)
1062

1063
        M = random_matrix(GF(2), k, n, sparse=False, density=p).row_space().basis_matrix()
1064
        M_n = random_matrix(GF(2), nwords, n, sparse=False, density=p)
1065
        B = LinearBinaryCodeStruct( M )
1066
        B_n = NonlinearBinaryCodeStruct( M_n )
1067
        B.run()
1068
        B_n.run()
1069

1070
        for i from 0 <= i < perms_per_code:
1071
            perm = [a-1 for a in list(S.random_element())]
1072
            C = LinearBinaryCodeStruct( matrix(GF(2), B.dimension, B.degree) )
1073
            C_n = NonlinearBinaryCodeStruct( matrix(GF(2), B_n.nwords, B_n.degree) )
1074
            for j from 0 <= j < B.dimension:
1075
                for h from 0 <= h < B.degree:
1076
                    bitset_set_to(&C.basis[j], perm[h], bitset_check(&B.basis[j], h))
1077
            for j from 0 <= j < B_n.nwords:
1078
                for h from 0 <= h < B_n.degree:
1079
                    bitset_set_to(&C_n.words[j], perm[h], bitset_check(&B_n.words[j], h))
1080
            # now C is a random permutation of B, and C_n of B_n
1081
            C.run()
1082
            C_n.run()
1083
            B_relab = B.canonical_relabeling()
1084
            C_relab = C.canonical_relabeling()
1085
            B_n_relab = B_n.canonical_relabeling()
1086
            C_n_relab = C_n.canonical_relabeling()
1087
            B_M = matrix(GF(2), B.dimension, B.degree)
1088
            C_M = matrix(GF(2), B.dimension, B.degree)
1089
            B_n_M = matrix(GF(2), B_n.nwords, B_n.degree)
1090
            C_n_M = matrix(GF(2), B_n.nwords, B_n.degree)
1091
            for j from 0 <= j < B.dimension:
1092
                for h from 0 <= h < B.degree:
1093
                    B_M[j,B_relab[h]] = bitset_check(&B.basis[j], h)
1094
                    C_M[j,C_relab[h]] = bitset_check(&C.basis[j], h)
1095
            for j from 0 <= j < B_n.nwords:
1096
                for h from 0 <= h < B_n.degree:
1097
                    B_n_M[j,B_n_relab[h]] = bitset_check(&B_n.words[j], h)
1098
                    C_n_M[j,C_n_relab[h]] = bitset_check(&C_n.words[j], h)
1099
            if B_M.row_space() != C_M.row_space():
1100
                print "can_lab error -- B:"
1101
                for j from 0 <= j < B.dimension:
1102
                    print bitset_string(&B.basis[j])
1103
                print perm
1104
                return
1105
            if sorted(B_n_M.rows()) != sorted(C_n_M.rows()):
1106
                print "can_lab error -- B_n:"
1107
                for j from 0 <= j < B_n.nwords:
1108
                    print bitset_string(&B_n.words[j])
1109
                print perm
1110
                return
1111
            isom = B.is_isomorphic(C)
1112
            if not isom:
1113
                print "isom -- B:"
1114
                for j from 0 <= j < B.dimension:
1115
                    print bitset_string(&B.basis[j])
1116
                print perm
1117
                print isom
1118
                return
1119
            isom = B_n.is_isomorphic(C_n)
1120
            if not isom:
1121
                print "isom -- B_n:"
1122
                for j from 0 <= j < B_n.nwords:
1123
                    print bitset_string(&B_n.words[j])
1124
                print perm
1125
                print isom
1126
                return
1127
                
1128
        num_tests += 4*perms_per_code
1129
        num_codes += 2
1130
        
1131
    print "All passed: %d random tests on %d codes."%(num_tests, num_codes)
1132

1133

1134

1135
Product

Resources

Company
