CoCalc -- refinement

GitHub Repository: sagemath/sagesmc
Path: blob/master/src/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
            sage_free(self.basis)
51
            sage_free(self.scratch_bitsets)
52
            sage_free(self.alpha_is_wd)
53
            PS_dealloc(self.word_ps)
54
            sage_free(self.alpha)
55
            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, NULL, 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
        output = <int *> sage_malloc(self.degree * sizeof(int))
310
        if part is NULL or ordering is NULL or output is NULL:
311
            PS_dealloc(part)
312
            sage_free(ordering)
313
            sage_free(output)
314
            raise MemoryError
315
        for i from 0 <= i < n:
316
            ordering[i] = i
317
        self.first_time = 1
318
        other.first_time = 1
319

320
        cdef bint isomorphic = double_coset(<void *> self, <void *> other, part, ordering, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_linear_codes, NULL, NULL, output)
321

322
        PS_dealloc(part)
323
        sage_free(ordering)
324
        if isomorphic:
325
            output_py = [output[i] for i from 0 <= i < n]
326
        else:
327
            output_py = False
328
        sage_free(output)
329
        return output_py
330

331
    def __dealloc__(self):
332
        cdef int j
333
        bitset_free(self.alpha_is_wd)
334
        for j from 0 <= j < 2*self.dimension+2:
335
            bitset_free(&self.scratch_bitsets[j])
336
        for j from 0 <= j < self.dimension:
337
            bitset_free(&self.basis[j])
338
        sage_free(self.basis); sage_free(self.scratch_bitsets)
339
        sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
340
        sage_free(self.alpha); sage_free(self.scratch)
341
        if self.output is not NULL:
342
            deallocate_agcl_output(self.output)
343

344
cdef int ith_word_linear(BinaryCodeStruct self, int i, bitset_s *word):
345
    cdef LinearBinaryCodeStruct LBCS = <LinearBinaryCodeStruct> self
346
    cdef int j
347
    bitset_zero(word)
348
    for j from 0 <= j < LBCS.dimension:
349
        if (1<<j)&i:
350
            bitset_xor(word, word, &LBCS.basis[j])
351
    return 0
352

353
cdef class NonlinearBinaryCodeStruct(BinaryCodeStruct):
354

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

366
        self.words = <bitset_s *> sage_malloc(self.nwords * sizeof(bitset_s))
367
        self.scratch_bitsets = <bitset_s *> sage_malloc((4*self.nwords+1) * sizeof(bitset_s))
368
        self.alpha_is_wd = <bitset_s *> sage_malloc(sizeof(bitset_s))
369
        self.word_ps = PS_new(self.nwords, 1)
370
        self.alpha = <int *> sage_malloc((self.nwords+self.degree) * sizeof(int))
371
        self.scratch = <int *> sage_malloc((3*self.nwords+3*self.degree+2) * sizeof(int))
372
        if self.words       is NULL or self.scratch_bitsets is NULL \
373
        or self.alpha_is_wd is NULL or self.word_ps         is NULL \
374
        or self.alpha       is NULL or self.scratch         is NULL:
375
            sage_free(self.words)
376
            sage_free(self.scratch_bitsets)
377
            sage_free(self.alpha_is_wd)
378
            PS_dealloc(self.word_ps)
379
            sage_free(self.alpha)
380
            sage_free(self.scratch)
381
            raise MemoryError
382

383
        cdef bint memerr = 0
384
        for i from 0 <= i < self.nwords:
385
            try: bitset_init(&self.words[i], self.degree)
386
            except MemoryError:
387
                for j from 0 <= j < i:
388
                    bitset_free(&self.words[j])
389
                memerr = 1
390
                break
391
        if not memerr:
392
            for i from 0 <= i < 4*self.nwords:
393
                try: bitset_init(&self.scratch_bitsets[i], self.degree)
394
                except MemoryError:
395
                    for j from 0 <= j < i:
396
                        bitset_free(&self.scratch_bitsets[j])
397
                    for j from 0 <= j < self.nwords:
398
                        bitset_free(&self.words[j])
399
                    memerr = 1
400
                    break
401
        if not memerr:
402
            try: bitset_init(&self.scratch_bitsets[4*self.nwords], self.nwords)
403
            except MemoryError:
404
                for j from 0 <= j < 4*self.nwords:
405
                    bitset_free(&self.scratch_bitsets[j])
406
                for j from 0 <= j < self.nwords:
407
                    bitset_free(&self.words[j])
408
                memerr = 1
409
        if not memerr:
410
            try: bitset_init(self.alpha_is_wd, self.nwords + self.degree)
411
            except MemoryError:
412
                for j from 0 <= j < 4*self.nwords + 1:
413
                    bitset_free(&self.scratch_bitsets[j])
414
                for j from 0 <= j < self.nwords:
415
                    bitset_free(&self.words[j])
416
                memerr = 1
417
        if memerr:
418
            sage_free(self.words); sage_free(self.scratch_bitsets)
419
            sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
420
            sage_free(self.alpha); sage_free(self.scratch)
421
            raise MemoryError
422
        else:
423
            bitset_zero(self.alpha_is_wd)
424
            for j from 0 <= j < self.nwords:
425
                bitset_zero(&self.words[j])
426

427
        if is_Matrix(arg):
428
            for i,j in arg.nonzero_positions():
429
                bitset_set(&self.words[i], j)
430

431
        self.output = NULL
432
        self.ith_word = &ith_word_nonlinear
433

434
    def __dealloc__(self):
435
        cdef int j
436
        bitset_free(self.alpha_is_wd)
437
        for j from 0 <= j < 4*self.nwords + 1:
438
            bitset_free(&self.scratch_bitsets[j])
439
        for j from 0 <= j < self.nwords:
440
            bitset_free(&self.words[j])
441
        sage_free(self.words); sage_free(self.scratch_bitsets)
442
        sage_free(self.alpha_is_wd); PS_dealloc(self.word_ps)
443
        sage_free(self.alpha); sage_free(self.scratch)
444
        if self.output is not NULL:
445
            deallocate_agcl_output(self.output)
446

447
    def run(self, partition=None):
448
        """
449
        Perform the canonical labeling and automorphism group computation,
450
        storing results to self.
451

452
        INPUT:
453
        partition -- an optional list of lists partition of the columns.
454
            default is the unit partition.
455

456
        EXAMPLES:
457
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
458

459
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,0,0,0],[0,0,1,0]]))
460
            sage: B.run()
461
            sage: B.automorphism_group()
462
            ([[2, 1, 0, 3], [0, 3, 2, 1]], 4, [1, 0])
463
            sage: B.canonical_relabeling()
464
            [2, 0, 3, 1]
465

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

473
            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]]))
474
            sage: B.run()
475
            sage: B.automorphism_group()
476
            ([[0, 1, 3, 2, 4, 5],
477
              [0, 2, 1, 3, 4, 5],
478
              [1, 0, 2, 3, 4, 5],
479
              [0, 1, 2, 3, 5, 4]],
480
             48,
481
             [4, 0, 1, 2])
482
            sage: B.canonical_relabeling()
483
            [2, 3, 4, 5, 0, 1]
484

485
        """
486
        cdef int n = self.degree
487
        cdef PartitionStack *part
488
        if partition is None:
489
            part = PS_new(n, 1)
490
        else:
491
            part = PS_from_list(partition)
492
        if part is NULL:
493
            raise MemoryError
494
        self.first_time = 1
495

496
        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, NULL, NULL)
497

498
        PS_dealloc(part)
499

500
    def automorphism_group(self):
501
        """
502
        Returns a list of generators of the automorphism group, along with its
503
        order and a base for which the list of generators is a strong generating
504
        set.
505

506
        EXAMPLE: (For more examples, see self.run())
507
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
508

509
            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]]))
510
            sage: B.run()
511
            sage: B.automorphism_group()
512
            ([[0, 1, 3, 2, 4, 5],
513
              [0, 2, 1, 3, 4, 5],
514
              [1, 0, 2, 3, 4, 5],
515
              [0, 1, 2, 3, 5, 4]],
516
             48,
517
             [4, 0, 1, 2])
518

519
        """
520
        cdef int i, j
521
        cdef list generators, base
522
        cdef Integer order
523
        if self.output is NULL:
524
            self.run()
525
        generators = []
526
        for i from 0 <= i < self.output.num_gens:
527
            generators.append([self.output.generators[i*self.degree + j] for j from 0 <= j < self.degree])
528
        order = Integer()
529
        SC_order(self.output.group, 0, order.value)
530
        base = [self.output.group.base_orbits[i][0] for i from 0 <= i < self.output.group.base_size]
531
        return generators, order, base
532

533
    def canonical_relabeling(self):
534
        """
535
        Returns a canonical relabeling (in list permutation format).
536

537
        EXAMPLES: (For more examples, see self.run())
538
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
539

540
            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]]))
541
            sage: B.run()
542
            sage: B.canonical_relabeling()
543
            [2, 3, 4, 5, 0, 1]
544

545
        """
546
        cdef int i
547
        if self.output is NULL:
548
            self.run()
549
        return [self.output.relabeling[i] for i from 0 <= i < self.degree]
550

551
    def is_isomorphic(self, NonlinearBinaryCodeStruct other):
552
        """
553
        Calculate whether self is isomorphic to other.
554

555
        EXAMPLES:
556
            sage: from sage.groups.perm_gps.partn_ref.refinement_binary import NonlinearBinaryCodeStruct
557

558
            sage: B = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]))
559
            sage: C = NonlinearBinaryCodeStruct(Matrix(GF(2), [[1,1,0,0,1,1],[1,1,1,1,0,0]]))
560
            sage: B.is_isomorphic(C)
561
            [2, 3, 0, 1, 4, 5]
562

563
        """
564
        cdef int i, n = self.degree
565
        cdef int *output, *ordering
566
        cdef PartitionStack *part
567
        part = PS_new(n, 1)
568
        ordering = <int *> sage_malloc(n * sizeof(int))
569
        output = <int *> sage_malloc(n * sizeof(int))
570
        if part is NULL or ordering is NULL or output is NULL:
571
            PS_dealloc(part)
572
            sage_free(ordering)
573
            sage_free(output)
574
            raise MemoryError
575
        for i from 0 <= i < n:
576
            ordering[i] = i
577
        self.first_time = 1
578
        other.first_time = 1
579

580
        cdef bint isomorphic = double_coset(<void *> self, <void *> other, part, ordering, n, &all_children_are_equivalent, &refine_by_bip_degree, &compare_nonlinear_codes, NULL, NULL, output)
581

582
        PS_dealloc(part)
583
        sage_free(ordering)
584
        if isomorphic:
585
            output_py = [output[i] for i from 0 <= i < n]
586
        else:
587
            output_py = False
588
        sage_free(output)
589
        return output_py
590

591
cdef int ith_word_nonlinear(BinaryCodeStruct self, int i, bitset_s *word):
592
    cdef NonlinearBinaryCodeStruct NBCS = <NonlinearBinaryCodeStruct> self
593
    bitset_copy(word, &NBCS.words[i])
594
    return 0
595

596
cdef int refine_by_bip_degree(PartitionStack *col_ps, void *S, int *cells_to_refine_by, int ctrb_len):
597
    """
598
    Refines the input partition by checking degrees of vertices to the given
599
    cells in the associated bipartite graph (vertices split into columns and
600
    words).
601

602
    INPUT:
603
    col_ps -- a partition stack, whose finest partition is the partition to be
604
        refined.
605
    S -- a binary code struct object
606
    cells_to_refine_by -- a list of pointers to cells to check degrees against
607
        in refining the other cells (updated in place)
608
    ctrb_len -- how many cells in cells_to_refine_by
609

610
    OUTPUT:
611

612
    An integer $I$ invariant under the orbits of $S_n$.  That is, if $\gamma$ is a
613
    permutation of the columns, then
614
    $$ I(G, PS, cells_to_refine_by) = I( \gamma(G), \gamma(PS), \gamma(cells_to_refine_by) ) .$$
615

616
    """
617
    cdef BinaryCodeStruct BCS = <BinaryCodeStruct> S
618
    cdef int current_cell_against = 0
619
    cdef int current_cell, i, r, j
620
    cdef int first_largest_subcell
621
    cdef int invariant = 0
622
    cdef PartitionStack *word_ps = BCS.word_ps
623
    cdef int *ctrb = BCS.alpha
624
    cdef bitset_s *ctrb_is_wd = BCS.alpha_is_wd
625

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

720
cdef int compare_linear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2, int degree):
721
    """
722
    Compare gamma_1(S1) and gamma_2(S2).
723

724
    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
725
    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
726

727
    Abstractly, what this function does is relabel the basis of B by gamma_1 and
728
    gamma_2, run a row reduction on each, and verify that the matrices are the
729
    same, which holds if and only if the rowspan is the same. In practice, if
730
    the codes are not equal, the reductions (which are carried out in an
731
    interleaved way) will terminate as soon as this is discovered, and whichever
732
    code has a 1 in the entry in which they differ is reported as larger.
733

734
    INPUT:
735
    gamma_1, gamma_2 -- list permutations (inverse)
736
    S1, S2 -- binary code struct objects
737

738
    """
739
    cdef int i, piv_loc_1, piv_loc_2, cur_col, cur_row=0
740
    cdef bint is_pivot_1, is_pivot_2
741
    cdef LinearBinaryCodeStruct BCS1 = <LinearBinaryCodeStruct> S1
742
    cdef LinearBinaryCodeStruct BCS2 = <LinearBinaryCodeStruct> S2
743
    cdef bitset_s *basis_1 = BCS1.scratch_bitsets                   # len = dim
744
    cdef bitset_s *basis_2 = &BCS1.scratch_bitsets[BCS1.dimension]   # len = dim
745
    cdef bitset_s *pivots  = &BCS1.scratch_bitsets[2*BCS1.dimension] # len 1
746
    cdef bitset_s *temp  = &BCS1.scratch_bitsets[2*BCS1.dimension+1] # len 1
747
    for i from 0 <= i < BCS1.dimension:
748
        bitset_copy(&basis_1[i], &BCS1.basis[i])
749
        bitset_copy(&basis_2[i], &BCS2.basis[i])
750
    bitset_zero(pivots)
751
    for cur_col from 0 <= cur_col < BCS1.degree:
752
        is_pivot_1 = 0
753
        is_pivot_2 = 0
754
        for i from cur_row <= i < BCS1.dimension:
755
            if bitset_check(&basis_1[i], gamma_1[cur_col]):
756
                is_pivot_1 = 1
757
                piv_loc_1 = i
758
                break
759
        for i from cur_row <= i < BCS1.dimension:
760
            if bitset_check(&basis_2[i], gamma_2[cur_col]):
761
                is_pivot_2 = 1
762
                piv_loc_2 = i
763
                break
764
        if is_pivot_1 != is_pivot_2:
765
            return <int>is_pivot_2 - <int>is_pivot_1
766
        if is_pivot_1:
767
            bitset_set(pivots, cur_col)
768
            if piv_loc_1 != cur_row:
769
                bitset_copy(temp, &basis_1[piv_loc_1])
770
                bitset_copy(&basis_1[piv_loc_1], &basis_1[cur_row])
771
                bitset_copy(&basis_1[cur_row], temp)
772
            if piv_loc_2 != cur_row:
773
                bitset_copy(temp, &basis_2[piv_loc_2])
774
                bitset_copy(&basis_2[piv_loc_2], &basis_2[cur_row])
775
                bitset_copy(&basis_2[cur_row], temp)
776
            for i from 0 <= i < cur_row:
777
                if bitset_check(&basis_1[i], gamma_1[cur_col]):
778
                    bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
779
                if bitset_check(&basis_2[i], gamma_2[cur_col]):
780
                    bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row])
781
            for i from cur_row < i < BCS1.dimension:
782
                if bitset_check(&basis_1[i], gamma_1[cur_col]):
783
                    bitset_xor(&basis_1[i], &basis_1[i], &basis_1[cur_row])
784
                if bitset_check(&basis_2[i], gamma_2[cur_col]):
785
                    bitset_xor(&basis_2[i], &basis_2[i], &basis_2[cur_row])
786
            cur_row += 1
787
        else:
788
            for i from 0 <= i < cur_row:
789
                if bitset_check(&basis_1[i], gamma_1[cur_col]) != bitset_check(&basis_2[i], gamma_2[cur_col]):
790
                    return <int>bitset_check(&basis_2[i], gamma_2[cur_col]) - <int>bitset_check(&basis_1[i], gamma_1[cur_col])
791
    return 0
792

793
cdef int compare_nonlinear_codes(int *gamma_1, int *gamma_2, void *S1, void *S2, int degree):
794
    """
795
    Compare gamma_1(S1) and gamma_2(S2).
796

797
    Return return -1 if gamma_1(S1) < gamma_2(S2), 0 if gamma_1(S1) == gamma_2(S2),
798
    1 if gamma_1(S1) > gamma_2(S2).  (Just like the python \code{cmp}) function.
799

800
    INPUT:
801
    gamma_1, gamma_2 -- list permutations (inverse)
802
    S1, S2 -- a binary code struct object
803

804
    """
805
    cdef int side=0, i, start, end, n_one_1, n_one_2, cur_col
806
    cdef int where_0, where_1
807
    cdef NonlinearBinaryCodeStruct BCS1 = <NonlinearBinaryCodeStruct> S1
808
    cdef NonlinearBinaryCodeStruct BCS2 = <NonlinearBinaryCodeStruct> S2
809
    cdef bitset_s *B_1_0 = BCS1.scratch_bitsets                   # nwords of len degree
810
    cdef bitset_s *B_1_1 = &BCS1.scratch_bitsets[BCS1.nwords]      # nwords of len degree
811
    cdef bitset_s *B_2_0 = &BCS1.scratch_bitsets[2*BCS1.nwords]    # nwords of len degree
812
    cdef bitset_s *B_2_1 = &BCS1.scratch_bitsets[3*BCS1.nwords]    # nwords of len degree
813
    cdef bitset_s *dividers = &BCS1.scratch_bitsets[4*BCS1.nwords] # 1 of len nwords
814
    cdef bitset_s *B_1_this, *B_1_other, *B_2_this, *B_2_other
815
    for i from 0 <= i < BCS1.nwords:
816
        bitset_copy(&B_1_0[i], &BCS1.words[i])
817
        bitset_copy(&B_2_0[i], &BCS2.words[i])
818
    bitset_zero(dividers)
819
    bitset_set(dividers, BCS1.nwords-1)
820

821
    for cur_col from 0 <= cur_col < BCS1.degree:
822
        if side == 0:
823
            B_1_this  = B_1_0
824
            B_1_other = B_1_1
825
            B_2_this  = B_2_0
826
            B_2_other = B_2_1
827
        else:
828
            B_1_this  = B_1_1
829
            B_1_other = B_1_0
830
            B_2_this  = B_2_1
831
            B_2_other = B_2_0
832
        side ^= 1
833
        start = 0
834
        while start < BCS1.nwords:
835
            end = start
836
            while not bitset_check(dividers, end):
837
                end += 1
838
            end += 1
839
            n_one_1 = 0
840
            n_one_2 = 0
841
            for i from start <= i < end:
842
                n_one_1 += bitset_check(&B_1_this[i], gamma_1[cur_col])
843
                n_one_2 += bitset_check(&B_2_this[i], gamma_2[cur_col])
844
            if n_one_1 != n_one_2:
845
                if n_one_1 > n_one_2:
846
                    return 1
847
                else:
848
                    return -1
849
            where_0 = start
850
            where_1 = end - n_one_1
851
            if start < where_1 and where_1 < end:
852
                bitset_set(dividers, where_1 - 1)
853
            for i from start <= i < end:
854
                if bitset_check(&B_1_this[i], gamma_1[cur_col]):
855
                    bitset_copy(&B_1_other[where_1], &B_1_this[i])
856
                    where_1 += 1
857
                else:
858
                    bitset_copy(&B_1_other[where_0], &B_1_this[i])
859
                    where_0 += 1
860
            where_0 = start
861
            where_1 = end - n_one_2
862
            for i from start <= i < end:
863
                if bitset_check(&B_2_this[i], gamma_2[cur_col]):
864
                    bitset_copy(&B_2_other[where_1], &B_2_this[i])
865
                    where_1 += 1
866
                else:
867
                    bitset_copy(&B_2_other[where_0], &B_2_this[i])
868
                    where_0 += 1
869
            start = end
870

871
    return 0
872

873
cdef bint all_children_are_equivalent(PartitionStack *col_ps, void *S):
874
    """
875
    Returns True if any refinement of the current partition results in the same
876
    structure.
877

878
    WARNING:
879
    Converse does not hold in general! See Lemma 2.25 of [1] for details, noting
880
    that the binary code is interpreted as a bipartite graph (see module docs
881
    for details).
882

883
    INPUT:
884
    col_ps -- the partition stack to be checked
885
    S -- a binary code struct object
886
    """
887
    cdef BinaryCodeStruct BCS = <BinaryCodeStruct> S
888
    cdef PartitionStack *word_ps = BCS.word_ps
889
    cdef int i, n = col_ps.degree + BCS.nwords
890
    cdef bint in_cell = 0
891
    cdef int nontrivial_cells = 0
892
    cdef int total_cells = PS_num_cells(col_ps) + PS_num_cells(word_ps)
893
    if n <= total_cells + 4:
894
        return 1
895
    for i from 0 <= i < BCS.nwords:
896
        if word_ps.levels[i] <= col_ps.depth:
897
            if in_cell:
898
                nontrivial_cells += 1
899
            in_cell = 0
900
        else:
901
            in_cell = 1
902
    in_cell = 0
903
    for i from 0 <= i < BCS.degree:
904
        if col_ps.levels[i] <= col_ps.depth:
905
            if in_cell:
906
                nontrivial_cells += 1
907
            in_cell = 0
908
        else:
909
            in_cell = 1
910
    if n == total_cells + nontrivial_cells:
911
        return 1
912
    if n == total_cells + nontrivial_cells + 1:
913
        return 1
914
    return 0
915

916
cdef inline int word_degree(PartitionStack *word_ps, BinaryCodeStruct BCS, int entry, int cell_index, PartitionStack *col_ps):
917
    """
918
    Returns the number of edges from the vertex corresponding to entry to
919
    vertices in the cell corresponding to cell_index.
920

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

946
cdef inline int col_degree(PartitionStack *col_ps, BinaryCodeStruct BCS, int entry, int cell_index, PartitionStack *word_ps):
947
    """
948
    Returns the number of edges from the vertex corresponding to entry to
949
    vertices in the cell corresponding to cell_index.
950

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

972
cdef inline int sort_by_function_codes(PartitionStack *PS, int start, int *degrees, int *counts, int *output, int count_max):
973
    """
974
    A simple counting sort, given the degrees of vertices to a certain cell.
975

976
    INPUT:
977
    PS -- the partition stack to be checked
978
    start -- beginning index of the cell to be sorted
979
    degrees -- the values to be sorted by
980
    count, count_max, output -- scratch space
981

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

1014

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

1020
    INPUT:
1021
    num -- run tests for this many codes
1022
    n_max -- test codes with at most this many columns
1023
    k_max -- test codes with at most this for dimension
1024
    perms_per_code -- test each code with this many random permutations
1025

1026
    DISCUSSION:
1027

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

1033
    For each code B (B2) generated, we uniformly generate perms_per_code random
1034
    permutations and verify that the canonical labels of B and the image of B
1035
    under the generated permutation are equal, and check that the double coset
1036
    function returns an isomorphism.
1037

1038
    TESTS::
1039

1040
        sage: import sage.groups.perm_gps.partn_ref.refinement_binary
1041
        sage: sage.groups.perm_gps.partn_ref.refinement_binary.random_tests()  # long time (up to 5s on sage.math, 2012)
1042
        All passed: ... random tests on ... codes.
1043

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

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

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

1125
        num_tests += 4*perms_per_code
1126
        num_codes += 2
1127

1128
    print "All passed: %d random tests on %d codes."%(num_tests, num_codes)
1129

1130

1131

1132
Product

Resources

Company
