Changeset 7760:39339387e58b


Ignore:
Timestamp:
11/11/07 12:23:36 (6 years ago)
Author:
R. L. Miller <rlmillster@…>
Branch:
default
Message:

3rd half

Location:
sage
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • sage/coding/binary_code.pxd

    r7758 r7760  
    33 
    44cdef class BinaryCode: 
    5     cdef unsigned int *columns 
    6     cdef int *ham_wts 
     5    cdef unsigned int *basis 
     6#    cdef unsigned int *columns 
     7    cdef unsigned int *words 
    78    cdef int ncols 
    89    cdef int nrows 
    … …  
    1314 
    1415cdef class OrbitPartition: 
     16    cdef unsigned int nwords 
     17    cdef int ncols 
    1518    cdef unsigned int *wd_parent 
    1619    cdef unsigned int *wd_rank 
    … …  
    2124    cdef int *col_min_cell_rep 
    2225    cdef int *col_size 
     26 
    2327    cdef unsigned int wd_find(self, unsigned int) 
    2428    cdef void wd_union(self, unsigned int, unsigned int) 
    2529    cdef int col_find(self, int) 
    2630    cdef void col_union(self, int, int) 
    27     cdef int merge_perm(self, int *, unsigned int *, int, int) 
     31    cdef int merge_perm(self, int *, unsigned int *) 
    2832 
    2933cdef class PartitionStack: 
    30     cdef int *wd_ents 
     34    cdef unsigned int *wd_ents 
    3135    cdef int *wd_lvls 
    3236    cdef int *col_ents 
    3337    cdef int *col_lvls 
    34     cdef int nwords 
     38    cdef unsigned int nwords 
    3539    cdef int ncols 
     40    cdef int radix 
     41    cdef unsigned int *col_degs   # 
     42    cdef int *col_counts          # 
     43    cdef int *col_output          # 
     44    cdef int *wd_degs             # 
     45    cdef unsigned int *wd_counts  # These are just for scratch space... 
     46    cdef unsigned int *wd_output  # 
     47 
    3648    cdef int is_discrete(self, int) 
    3749    cdef int num_cells(self, int) 
    3850    cdef int sat_225(self, int) 
    39     cdef int is_min_cell_rep(self, int, int, int) 
    40     cdef int is_fixed(self, int, int, int) 
    41     cdef void col_percolate(self, int start, int end) 
    42     cdef void wd_percolate(self, int start, int end) 
    43     cdef int split_vertex(self, int, int, int) 
    44     cdef int col_degree(self, BinaryCode, int, int, int) 
    45     cdef int wd_degree(self, BinaryCode, int, int, int) 
    46     cdef int sort_cols(self, int, int *, int) 
    47     cdef int sort_wds(self, int, int *, int) 
    48     cdef int refine(self, int, int *, int *, BinaryCode) 
     51    cdef unsigned int min_cell_reps(self, int) 
     52    cdef unsigned int fixed_cols(self, unsigned int, int) 
     53    cdef unsigned int first_smallest_nontrivial(self, int) 
     54    cdef void col_percolate(self, int, int) 
     55    cdef void wd_percolate(self, unsigned int, unsigned int) 
     56    cdef int split_column(self, int, int) 
     57    cdef unsigned int col_degree(self, BinaryCode, int, unsigned int, int) 
     58    cdef int wd_degree(self, BinaryCode, unsigned int, int, int) 
     59    cdef int sort_cols(self, int, int) 
     60    cdef unsigned int sort_wds(self, unsigned int, int) 
     61 
     62################################################################################ 
     63################################################################################ 
     64################################################################################ 
     65 
     66    cdef unsigned int refine(self, int, unsigned int *, int *, BinaryCode) 
    4967    cdef void get_permutation(self, PartitionStack, int *, int *) 
    5068    cdef int cmp(self, PartitionStack, BinaryCode) 
     69 
     70cdef class BinaryCodeClassifier: 
     71    cdef int *ham_wts 
     72 
     73 
     74 
     75 
     76 
     77 
     78 
     79 
     80 
     81 
     82 
     83 
     84 
     85 
     86 
     87 
     88 
     89 

  • sage/coding/binary_code.pyx

    r7758 r7760  
    1717AUTHOR: 
    1818    Robert L Miller (Oct-Nov 2007) 
    19         * Compiled code datastructure 
     19        * compiled code datastructure 
    2020        * union-find based orbit partition 
    2121        * optimized partition stack class 
    2222 
     23        * NICE-based partition refinement algorithm 
     24        * canonical generation function 
     25 
    2326""" 
    2427 
    25 #***************************************************************************** 
     28#******************************************************************************* 
    2629#         Copyright (C) 2007 Robert L. Miller <rlmillster@gmail.com> 
    2730# 
    2831# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL) 
    2932#                         http://www.gnu.org/licenses/ 
    30 #***************************************************************************** 
     33#******************************************************************************* 
    3134 
    3235include '../ext/cdefs.pxi' 
    … …  
    3841from sage.rings.integer import Integer 
    3942 
     43## NOTE - Since most of the functions are used from within the module, cdef'd 
     44## functions come without an underscore, and the def'd equivalents, which are 
     45## essentially only for doctesting and debugging, have underscores. 
     46 
    4047cdef class BinaryCode: 
    4148    """ 
    4249    Minimal, but optimized, binary code object. 
    4350 
    44     EXAMPLES: 
    45         sage: import sage.coding.binary_code 
    46         sage: from sage.coding.binary_code import * 
    47         sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]]) 
    48         sage: B = BinaryCode(M) 
    49         sage: B 
    50         Binary [8,4] linear code, generator matrix 
    51         [11110000] 
    52         [00111100] 
    53         [00001111] 
    54         [10101010] 
    55         sage: B._is_one(7, 4) 
    56         0 
    57         sage: B._is_one(8, 4) 
    58         1 
    59         sage: B._is_automorphism([1,0,3,2,4,5,6,7], [1, 2, 4, 9]) 
    60         1 
    61  
    6251    """ 
    6352    def __new__(self, arg1): 
    64         """ 
    65         Create binary code from matrix. Input is assumed to be a matrix 
    66         over $GF(2)$. 
    67  
    68         EXAMPLE: 
    69         sage: import sage.coding.binary_code 
    70         sage: from sage.coding.binary_code import * 
    71         sage: M = Matrix(GF(2), [[1,1,1,1,0,0,0,0],[0,0,1,1,1,1,0,0],[0,0,0,0,1,1,1,1],[1,0,1,0,1,0,1,0]]) 
    72         sage: B = BinaryCode(M) 
    73         sage: B 
    74         Binary [8,4] linear code, generator matrix 
    75         [11110000] 
    76         [00111100] 
    77         [00001111] 
    78         [10101010] 
    79  
    80         """ 
    81         cdef int i, j 
    82         cdef unsigned int k 
     53        cdef unsigned int parity, word, combination 
     54        cdef int nrows, i, j 
     55        cdef unsigned int *self_words, *self_basis 
    8356        self.radix = 8*sizeof(unsigned int) 
    8457        self.ncols = arg1.ncols() 
    … …  
    8659        self.nwords = 1 << self.nrows 
    8760 
    88         if self.nrows > self.radix or ceil(log(self.ncols,2)) > self.radix: 
     61        if self.nrows > self.radix or self.ncols > self.radix: 
    8962            raise NotImplementedError("Columns and rows are stored as ints. This code is too big.") 
    90         self.columns = <unsigned int *> sage_malloc( self.ncols * sizeof(unsigned int) ) 
    91         if not self.columns: 
     63 
     64        self.words = <unsigned int *> sage_malloc( self.nwords * sizeof(unsigned int) ) 
     65        self.basis = <unsigned int *> sage_malloc( self.nrows * sizeof(unsigned int) ) 
     66        if not self.words or not self.basis: 
     67            if self.words: sage_free(self.words) 
     68            if self.basis: sage_free(self.basis) 
    9269            raise MemoryError("Memory.") 
    9370 
    94         cols = arg1.columns() 
    95         for i from 0 <= i < self.ncols: 
    96             k = 0 
    97             for j in cols[i].nonzero_positions(): 
    98                 k += (1 << j) 
    99             self.columns[i] = k 
    100  
    101         self.ham_wts = <int *> sage_malloc( 16 * sizeof(int) ) 
    102         if not self.ham_wts: 
    103             sage_free(self.columns) 
    104             raise MemoryError("Memory.") 
    105         self.ham_wts[0]  = 0; self.ham_wts[1]  = 1; self.ham_wts[2]  = 1; self.ham_wts[3]  = 2 
    106         self.ham_wts[4]  = 1; self.ham_wts[5]  = 2; self.ham_wts[6]  = 2; self.ham_wts[7]  = 3 
    107         self.ham_wts[8]  = 1; self.ham_wts[9]  = 2; self.ham_wts[10] = 2; self.ham_wts[11] = 3 
    108         self.ham_wts[12] = 2; self.ham_wts[13] = 3; self.ham_wts[14] = 3; self.ham_wts[15] = 4 
     71        nrows = self.nrows 
     72        self_words = self.words 
     73        self_basis = self.basis 
     74 
     75        rows = arg1.rows() 
     76        for i from 0 <= i < nrows: 
     77            word = 0 
     78            for j in rows[i].nonzero_positions(): 
     79                word += (1<<j) 
     80            self_basis[i] = word 
     81 
     82        word = 0 
     83        parity = 0 
     84        combination = 0 
     85        while True: 
     86            self_words[combination] = word 
     87            parity ^= 1 
     88            j = 0 
     89            if not parity: 
     90                while not combination & (1 << j): j += 1 
     91                j += 1 
     92            if j == nrows: break 
     93            else: 
     94                combination ^= (1 << j) 
     95                word ^= self_basis[j] 
     96        # NOTE: when implementing construction from parent 
     97        # matrices, remember to grab this data from the parent 
    10998 
    11099    def __dealloc__(self): 
    111         sage_free(self.columns) 
    112         sage_free(self.ham_wts) 
     100        sage_free(self.words) 
     101        sage_free(self.basis) 
     102 
     103    def delete(self): 
     104        sage_free(self.words) 
     105        sage_free(self.basis) 
    113106 
    114107    def __repr__(self): 
    … …  
    166159 
    167160    cdef int is_one(self, unsigned int word, int column): 
    168         cdef int i, j 
    169         cdef unsigned int k 
    170         i = 0 
    171         k = word & self.columns[column] 
    172         for j from 0 <= j < self.radix by 4: 
    173             i += self.ham_wts[15 & k] 
    174             k = k >> 4 
    175         return i % 2 
     161        if self.words[word] & (1 << column): 
     162            return 1 
     163        else: 
     164            return 0 
    176165 
    177166    def _is_automorphism(self, col_gamma, basis_gamma): 
    … …  
    221210 
    222211    cdef int is_automorphism(self, int *col_gamma, unsigned int *basis_gamma): 
    223         cdef int i, j 
    224         for i from 0 <= i < self.nrows: 
    225             for j from 0 <= j < self.ncols: 
     212        cdef int i, j, self_nrows = self.nrows, self_ncols = self.ncols 
     213        for i from 0 <= i < self_nrows: 
     214            for j from 0 <= j < self_ncols: 
    226215                if self.is_one(1 << i, j) != self.is_one(basis_gamma[i], col_gamma[j]): 
    227216                    return 0 
    … …  
    243232        cdef int col 
    244233        nwords = (1 << nrows) 
     234        self.nwords = nwords 
     235        self.ncols = ncols 
    245236        self.wd_parent =       <unsigned int *> sage_malloc( nwords * sizeof(unsigned int) ) 
    246237        self.wd_rank =         <unsigned int *> sage_malloc( nwords * sizeof(unsigned int) ) 
    … …  
    282273        sage_free(self.col_size) 
    283274 
     275    def delete(self): 
     276        sage_free(self.wd_parent) 
     277        sage_free(self.wd_rank) 
     278        sage_free(self.wd_min_cell_rep) 
     279        sage_free(self.wd_size) 
     280        sage_free(self.col_parent) 
     281        sage_free(self.col_rank) 
     282        sage_free(self.col_min_cell_rep) 
     283        sage_free(self.col_size) 
     284 
     285    def __repr__(self): 
     286        """ 
     287        Return a string representation of the orbit partition. 
     288 
     289        EXAMPLE: 
     290            sage: import sage.coding.binary_code 
     291            sage: from sage.coding.binary_code import * 
     292            sage: O = OrbitPartition(4, 8) 
     293            sage: O 
     294            OrbitPartition on 16 words and 8 columns. Data: 
     295            Words: 
     296            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     297            Columns: 
     298            0,1,2,3,4,5,6,7 
     299 
     300        """ 
     301        cdef unsigned int i 
     302        cdef int j 
     303        s = 'OrbitPartition on %d words and %d columns. Data:\nWords:\n'%(self.nwords, self.ncols) 
     304        for i from 0 <= i < self.nwords: 
     305            s += '%d,'%self.wd_parent[i] 
     306        s = s[:-1] + '\nColumns:\n' 
     307        for j from 0 <= j < self.ncols: 
     308            s += '%d,'%self.col_parent[j] 
     309        return s[:-1] 
     310 
     311    def _wd_find(self, word): 
     312        """ 
     313        Returns the root of word. 
     314 
     315        EXAMPLE: 
     316            sage: import sage.coding.binary_code 
     317            sage: from sage.coding.binary_code import * 
     318            sage: O = OrbitPartition(4, 8) 
     319            sage: O 
     320            OrbitPartition on 16 words and 8 columns. Data: 
     321            Words: 
     322            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     323            Columns: 
     324            0,1,2,3,4,5,6,7 
     325            sage: O._wd_find(12) 
     326            12L 
     327 
     328        """ 
     329        return self.wd_find(word) 
     330 
    284331    cdef unsigned int wd_find(self, unsigned int word): 
    285332        if self.wd_parent[word] == word: 
    … …  
    288335            self.wd_parent[word] = self.wd_find(self.wd_parent[word]) 
    289336            return self.wd_parent[word] 
     337 
     338    def _wd_union(self, x, y): 
     339        """ 
     340        Join the cells containing x and y. 
     341 
     342        EXAMPLE: 
     343            sage: import sage.coding.binary_code 
     344            sage: from sage.coding.binary_code import * 
     345            sage: O = OrbitPartition(4, 8) 
     346            sage: O 
     347            OrbitPartition on 16 words and 8 columns. Data: 
     348            Words: 
     349            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     350            Columns: 
     351            0,1,2,3,4,5,6,7 
     352            sage: O._wd_union(1,10) 
     353            sage: O 
     354            OrbitPartition on 16 words and 8 columns. Data: 
     355            Words: 
     356            0,1,2,3,4,5,6,7,8,9,1,11,12,13,14,15 
     357            Columns: 
     358            0,1,2,3,4,5,6,7 
     359            sage: O._wd_find(10) 
     360            1L 
     361 
     362        """ 
     363        self.wd_union(x, y) 
    290364 
    291365    cdef void wd_union(self, unsigned int x, unsigned int y): 
    … …  
    307381            self.wd_rank[x_root] += 1 
    308382 
     383    def _col_find(self, col): 
     384        """ 
     385        Returns the root of col. 
     386 
     387        EXAMPLE: 
     388            sage: import sage.coding.binary_code 
     389            sage: from sage.coding.binary_code import * 
     390            sage: O = OrbitPartition(4, 8) 
     391            sage: O 
     392            OrbitPartition on 16 words and 8 columns. Data: 
     393            Words: 
     394            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     395            Columns: 
     396            0,1,2,3,4,5,6,7 
     397            sage: O._col_find(6) 
     398            6 
     399 
     400        """ 
     401        return self.col_find(col) 
     402 
    309403    cdef int col_find(self, int col): 
    310404        if self.col_parent[col] == col: 
    … …  
    313407            self.col_parent[col] = self.col_find(self.col_parent[col]) 
    314408            return self.col_parent[col] 
     409 
     410    def _col_union(self, x, y): 
     411        """ 
     412        Join the cells containing x and y. 
     413 
     414        EXAMPLE: 
     415            sage: import sage.coding.binary_code 
     416            sage: from sage.coding.binary_code import * 
     417            sage: O = OrbitPartition(4, 8) 
     418            sage: O 
     419            OrbitPartition on 16 words and 8 columns. Data: 
     420            Words: 
     421            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     422            Columns: 
     423            0,1,2,3,4,5,6,7 
     424            sage: O._col_union(1,4) 
     425            sage: O 
     426            OrbitPartition on 16 words and 8 columns. Data: 
     427            Words: 
     428            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     429            Columns: 
     430            0,1,2,3,1,5,6,7 
     431            sage: O._col_find(4) 
     432            1 
     433 
     434        """ 
     435        self.col_union(x, y) 
    315436 
    316437    cdef void col_union(self, int x, int y): 
    … …  
    332453            self.col_rank[x_root] += 1 
    333454 
    334     cdef int merge_perm(self, int *col_gamma, unsigned int *wd_gamma, int nrows, int ncols): 
    335         """ 
    336         Returns 0 if nothing was done, otherwise returns 1. 
    337  
    338         If gamma[a] = b, then after merge_perm, a and b will be in the same cell. 
    339  
    340         """ 
     455    def _merge_perm(self, col_gamma, wd_gamma): 
     456        """ 
     457        Merges the cells of self under the given permutation. If gamma[a] = b, 
     458        then after merge_perm, a and b will be in the same cell. Returns 0 if 
     459        nothing was done, otherwise returns 1. 
     460 
     461        EXAMPLE: 
     462            sage: import sage.coding.binary_code 
     463            sage: from sage.coding.binary_code import * 
     464            sage: O = OrbitPartition(4, 8) 
     465            sage: O 
     466            OrbitPartition on 16 words and 8 columns. Data: 
     467            Words: 
     468            0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
     469            Columns: 
     470            0,1,2,3,4,5,6,7 
     471            sage: O._merge_perm([1,0,3,2,4,5,6,7], [0,1,2,3,4,5,6,7,9,8,11,10,13,12,15,14]) 
     472            1 
     473            sage: O 
     474            OrbitPartition on 16 words and 8 columns. Data: 
     475            Words: 
     476            0,1,2,3,4,5,6,7,8,8,10,10,12,12,14,14 
     477            Columns: 
     478            0,0,2,2,4,5,6,7 
     479 
     480        """ 
     481        cdef unsigned int i 
     482        cdef int *_col_gamma 
     483        cdef unsigned int *_wd_gamma 
     484        _wd_gamma = <unsigned int *> sage_malloc(self.nwords * sizeof(unsigned int)) 
     485        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int)) 
     486        if not (_col_gamma and _wd_gamma): 
     487            if _wd_gamma: sage_free(_wd_gamma) 
     488            if _col_gamma: sage_free(_col_gamma) 
     489            raise MemoryError("Memory.") 
     490        for i from 0 <= i < self.nwords: 
     491            _wd_gamma[i] = wd_gamma[i] 
     492        for i from 0 <= i < self.ncols: 
     493            _col_gamma[i] = col_gamma[i] 
     494        result = self.merge_perm(_col_gamma, _wd_gamma) 
     495        sage_free(_col_gamma) 
     496        sage_free(_wd_gamma) 
     497        return result 
     498 
     499    cdef int merge_perm(self, int *col_gamma, unsigned int *wd_gamma): 
    341500        cdef unsigned int i, gamma_i_root 
    342501        cdef int j, gamma_j_root, return_value = 0 
    343         for i from 0 <= i < (1 << nrows): 
    344             if self.wd_parent[i] == i: 
     502        cdef unsigned int *self_wd_parent = self.wd_parent 
     503        cdef int *self_col_parent = self.col_parent 
     504        for i from 0 <= i < self.nwords: 
     505            if self_wd_parent[i] == i: 
    345506                gamma_i_root = self.wd_find(wd_gamma[i]) 
    346507                if gamma_i_root != i: 
    347508                    return_value = 1 
    348509                    self.wd_union(i, gamma_i_root) 
    349         for j from 0 <= j < ncols: 
    350             if self.col_parent[j] == j: 
     510        for j from 0 <= j < self.ncols: 
     511            if self_col_parent[j] == j: 
    351512                gamma_j_root = self.col_find(col_gamma[j]) 
    352513                if gamma_j_root != j: 
    … …  
    356517 
    357518cdef class PartitionStack: 
    358  
     519    """ 
     520    Partition stack structure for traversing the search tree during automorphism 
     521    group computation. 
     522 
     523    """ 
    359524    def __new__(self, arg1, arg2=None): 
    360525        cdef int k 
    361526        self.nwords = int(arg1) 
    362527        self.ncols = int(arg2) 
    363         self.wd_ents = <int *> sage_malloc( self.nwords * sizeof(int) ) 
    364         self.wd_lvls = <int *> sage_malloc( self.nwords * sizeof(int) ) 
    365         self.col_ents = <int *> sage_malloc( self.ncols * sizeof(int) ) 
    366         self.col_lvls = <int *> sage_malloc( self.ncols * sizeof(int) ) 
    367         if not (self.wd_ents and self.wd_lvls and self.col_ents and self.col_lvls): 
     528        self.radix      = 8*sizeof(unsigned int) 
     529 
     530        # data 
     531        self.wd_ents    = <unsigned int *> sage_malloc( self.nwords * sizeof(unsigned int) ) 
     532        self.wd_lvls    = <int *> sage_malloc( self.nwords * sizeof(int) ) 
     533        self.col_ents   = <int *> sage_malloc( self.ncols * sizeof(int) ) 
     534        self.col_lvls   = <int *> sage_malloc( self.ncols * sizeof(int) ) 
     535 
     536        # scratch space 
     537        self.col_degs   = <unsigned int *> sage_malloc( self.ncols * sizeof(unsigned int) ) 
     538        self.col_counts = <int *> sage_malloc( self.nwords * sizeof(int) ) 
     539        self.col_output = <int *> sage_malloc( self.ncols * sizeof(int) ) 
     540        self.wd_degs    = <int *> sage_malloc( self.nwords * sizeof(int) ) 
     541        self.wd_counts  = <unsigned int *> sage_malloc( self.ncols * sizeof(unsigned int) ) 
     542        self.wd_output  = <unsigned int *> sage_malloc( self.nwords * sizeof(unsigned int) ) 
     543 
     544        if not (self.wd_ents and self.wd_lvls and self.col_ents and self.col_lvls \ 
     545            and self.col_degs and self.col_counts and self.col_output \ 
     546            and self.wd_degs and self.wd_counts and self.wd_output): 
    368547            if self.wd_ents: sage_free(self.wd_ents) 
    369548            if self.wd_lvls: sage_free(self.wd_lvls) 
    370549            if self.col_ents: sage_free(self.col_ents) 
    371550            if self.col_lvls: sage_free(self.col_lvls) 
     551            if self.col_degs: sage_free(self.col_degs) 
     552            if self.col_counts: sage_free(self.col_counts) 
     553            if self.col_output: sage_free(self.col_output) 
     554            if self.wd_degs: sage_free(self.wd_degs) 
     555            if self.wd_counts: sage_free(self.wd_counts) 
     556            if self.wd_output: sage_free(self.wd_output) 
    372557            raise MemoryError("Memory.") 
     558 
    373559        for k from 0 <= k < self.nwords-1: 
    374560            self.wd_ents[k] = k 
    … …  
    387573        sage_free(self.col_ents) 
    388574        sage_free(self.col_lvls) 
     575        sage_free(self.col_degs) 
     576        sage_free(self.col_counts) 
     577        sage_free(self.col_output) 
     578        sage_free(self.wd_degs) 
     579        sage_free(self.wd_counts) 
     580        sage_free(self.wd_output) 
     581 
     582    def delete(self): 
     583        sage_free(self.wd_ents) 
     584        sage_free(self.wd_lvls) 
     585        sage_free(self.col_ents) 
     586        sage_free(self.col_lvls) 
     587        sage_free(self.col_degs) 
     588        sage_free(self.col_counts) 
     589        sage_free(self.col_output) 
     590        sage_free(self.wd_degs) 
     591        sage_free(self.wd_counts) 
     592        sage_free(self.wd_output) 
    389593 
    390594    def __repr__(self): 
     595        """ 
     596        Return a string representation of self. 
     597 
     598        EXAMPLE: 
     599            sage: import sage.coding.binary_code 
     600            sage: from sage.coding.binary_code import * 
     601            sage: P = PartitionStack(4, 6) 
     602            sage: P 
     603            ({0,1,2,3}) 
     604            ({0,1,2,3}) 
     605            ({0,1,2,3}) 
     606            ({0,1,2,3}) 
     607            ({0,1,2,3,4,5}) 
     608            ({0,1,2,3,4,5}) 
     609            ({0,1,2,3,4,5}) 
     610            ({0,1,2,3,4,5}) 
     611            ({0,1,2,3,4,5}) 
     612            ({0,1,2,3,4,5}) 
     613 
     614        """ 
    391615        cdef int i, j, k 
    392616        s = '' 
    … …  
    411635        return s 
    412636 
     637    def _is_discrete(self, k): 
     638        """ 
     639        Returns whether the partition at level k is discrete. 
     640         
     641        EXAMPLE: 
     642            sage: import sage.coding.binary_code 
     643            sage: from sage.coding.binary_code import * 
     644            sage: P = PartitionStack(4, 6) 
     645            sage: [P._split_column(i,i+1) for i in range(5)] 
     646            [0, 1, 2, 3, 4] 
     647            sage: P 
     648            ({0,1,2,3}) 
     649            ({0,1,2,3}) 
     650            ({0,1,2,3}) 
     651            ({0,1,2,3}) 
     652            ({0,1,2,3,4,5}) 
     653            ({0},{1,2,3,4,5}) 
     654            ({0},{1},{2,3,4,5}) 
     655            ({0},{1},{2},{3,4,5}) 
     656            ({0},{1},{2},{3},{4,5}) 
     657            ({0},{1},{2},{3},{4},{5}) 
     658            sage: P._is_discrete(4) 
     659            0 
     660            sage: P._is_discrete(5) 
     661            1 
     662 
     663        """ 
     664        return self.is_discrete(k) 
     665 
    413666    cdef int is_discrete(self, int k): 
    414667        cdef int i 
     668        cdef int *self_col_lvls = self.col_lvls 
    415669        for i from 0 <= i < self.ncols: 
    416             if self.col_lvls[i] > k: 
     670            if self_col_lvls[i] > k: 
    417671                return 0 
    418672        return 1 
    419673 
     674    def _num_cells(self, k): 
     675        """ 
     676        Returns the number of cells in the partition at level k. 
     677 
     678        EXAMPLE: 
     679            sage: import sage.coding.binary_code 
     680            sage: from sage.coding.binary_code import * 
     681            sage: P = PartitionStack(4, 6) 
     682            sage: [P._split_column(i,i+1) for i in range(5)] 
     683            [0, 1, 2, 3, 4] 
     684            sage: P 
     685            ({0,1,2,3}) 
     686            ({0,1,2,3}) 
     687            ({0,1,2,3}) 
     688            ({0,1,2,3}) 
     689            ({0,1,2,3,4,5}) 
     690            ({0},{1,2,3,4,5}) 
     691            ({0},{1},{2,3,4,5}) 
     692            ({0},{1},{2},{3,4,5}) 
     693            ({0},{1},{2},{3},{4,5}) 
     694            ({0},{1},{2},{3},{4},{5}) 
     695            sage: P._num_cells(3) 
     696            5 
     697 
     698        """ 
     699        return self.num_cells(k) 
     700 
    420701    cdef int num_cells(self, int k): 
    421702        cdef int i, j = 0 
     703        cdef int *self_wd_lvls = self.wd_lvls 
     704        cdef int *self_col_lvls = self.col_lvls 
    422705        for i from 0 <= i < self.nwords: 
    423             if self.wd_lvls[i] <= k: 
     706            if self_wd_lvls[i] <= k: 
    424707                j += 1 
    425708        for i from 0 <= i < self.ncols: 
    426             if self.col_lvls[i] <= k: 
     709            if self_col_lvls[i] <= k: 
    427710                j += 1 
    428711        return j 
     712 
     713    def _sat_225(self, k): 
     714        """ 
     715        Returns whether the partition at level k satisfies the hypotheses of 
     716        Lemma 2.25 in Brendan McKay's Practical Graph Isomorphism paper (see 
     717        sage/graphs/graph_isom.pyx. 
     718 
     719        EXAMPLE: 
     720            sage: import sage.coding.binary_code 
     721            sage: from sage.coding.binary_code import * 
     722            sage: P = PartitionStack(4, 6) 
     723            sage: [P._split_column(i,i+1) for i in range(5)] 
     724            [0, 1, 2, 3, 4] 
     725            sage: P._sat_225(3) 
     726            0 
     727            sage: P._sat_225(4) 
     728            1 
     729            sage: P 
     730            ({0,1,2,3}) 
     731            ({0,1,2,3}) 
     732            ({0,1,2,3}) 
     733            ({0,1,2,3}) 
     734            ({0,1,2,3,4,5}) 
     735            ({0},{1,2,3,4,5}) 
     736            ({0},{1},{2,3,4,5}) 
     737            ({0},{1},{2},{3,4,5}) 
     738            ({0},{1},{2},{3},{4,5}) 
     739            ({0},{1},{2},{3},{4},{5}) 
     740 
     741        """ 
     742        return self.sat_225(k) 
    429743 
    430744    cdef int sat_225(self, int k): 
    431745        cdef int i, n = self.nwords + self.ncols, in_cell = 0 
    432746        cdef int nontrivial_cells = 0, total_cells = self.num_cells(k) 
     747        cdef int *self_wd_lvls = self.wd_lvls 
     748        cdef int *self_col_lvls = self.col_lvls 
    433749        if n <= total_cells + 4: 
    434750            return 1 
    435751        for i from 0 <= i < self.nwords: 
    436             if self.wd_lvls[i] <= k: 
     752            if self_wd_lvls[i] <= k: 
    437753                if in_cell: 
    438754                    nontrivial_cells += 1 
    … …  
    442758        in_cell = 0 
    443759        for i from 0 <= i < self.ncols: 
    444             if self.col_lvls[i] <= k: 
     760            if self_col_lvls[i] <= k: 
    445761                if in_cell: 
    446762                    nontrivial_cells += 1 
    … …  
    454770        return 0 
    455771 
    456     cdef int is_min_cell_rep(self, int col_not_wd, int i, int k): 
    457         if i == 0: 
    458             return 1 
    459         if col_not_wd: 
    460             return self.col_lvls[i-1] <= k 
    461         else: 
    462             return self.wd_lvls[i-1] <= k 
    463  
    464     cdef int is_fixed(self, int col_not_wd, int i, int k): 
    465         """ 
    466         Assuming it is a min cell rep. 
    467         """ 
    468         if col_not_wd: 
    469             return self.col_lvls[i] <= k 
    470         else: 
    471             return self.wd_lvls[i] <= k 
    472  
    473     ## TODO: first smallest nontrivial 
    474  
     772    def _min_cell_reps(self, k): 
     773        """ 
     774        Returns an integer whose bits represent which columns are minimal cell 
     775        representatives. 
     776 
     777        EXAMPLE: 
     778            sage: import sage.coding.binary_code 
     779            sage: from sage.coding.binary_code import * 
     780            sage: P = PartitionStack(4, 6) 
     781            sage: [P._split_column(i,i+1) for i in range(5)] 
     782            [0, 1, 2, 3, 4] 
     783            sage: a = P._min_cell_reps(2) 
     784            sage: Integer(a).binary() 
     785            '111' 
     786            sage: P 
     787            ({0,1,2,3}) 
     788            ({0,1,2,3}) 
     789            ({0,1,2,3}) 
     790            ({0,1,2,3}) 
     791            ({0,1,2,3,4,5}) 
     792            ({0},{1,2,3,4,5}) 
     793            ({0},{1},{2,3,4,5}) 
     794            ({0},{1},{2},{3,4,5}) 
     795            ({0},{1},{2},{3},{4,5}) 
     796            ({0},{1},{2},{3},{4},{5}) 
     797 
     798        """ 
     799        return self.min_cell_reps(k) 
     800 
     801    cdef unsigned int min_cell_reps(self, int k): 
     802        cdef int i 
     803        cdef unsigned int reps = 1 
     804        cdef int *self_col_lvls = self.col_lvls 
     805        for i from 0 < i < self.ncols: 
     806            if self_col_lvls[i-1] <= k: 
     807                reps += (1 << i) 
     808        return reps 
     809 
     810    def _fixed_cols(self, mcrs, k): 
     811        """ 
     812        Returns an integer whose bits represent which columns are fixed. For 
     813        efficiency, mcrs is the output of min_cell_reps. 
     814 
     815        EXAMPLE: 
     816            sage: import sage.coding.binary_code 
     817            sage: from sage.coding.binary_code import * 
     818            sage: P = PartitionStack(4, 6) 
     819            sage: [P._split_column(i,i+1) for i in range(5)] 
     820            [0, 1, 2, 3, 4] 
     821            sage: a = P._fixed_cols(7, 2) 
     822            sage: Integer(a).binary() 
     823            '11' 
     824            sage: P 
     825            ({0,1,2,3}) 
     826            ({0,1,2,3}) 
     827            ({0,1,2,3}) 
     828            ({0,1,2,3}) 
     829            ({0,1,2,3,4,5}) 
     830            ({0},{1,2,3,4,5}) 
     831            ({0},{1},{2,3,4,5}) 
     832            ({0},{1},{2},{3,4,5}) 
     833            ({0},{1},{2},{3},{4,5}) 
     834            ({0},{1},{2},{3},{4},{5}) 
     835 
     836        """ 
     837        return self.fixed_cols(mcrs, k) 
     838 
     839    cdef unsigned int fixed_cols(self, unsigned int mcrs, int k): 
     840        cdef int i 
     841        cdef unsigned int fixed = 0 
     842        cdef int *self_col_lvls = self.col_lvls 
     843        for i from 0 <= i < self.ncols: 
     844            if self_col_lvls[i] <= k: 
     845                fixed += (1 << i) 
     846        return fixed & mcrs 
     847 
     848    def _first_smallest_nontrivial(self, k): 
     849        """ 
     850        Returns an integer representing the first, smallest nontrivial cell. 
     851 
     852        EXAMPLE: 
     853            sage: import sage.coding.binary_code 
     854            sage: from sage.coding.binary_code import * 
     855            sage: P = PartitionStack(4, 6) 
     856            sage: [P._split_column(i,i+1) for i in range(5)] 
     857            [0, 1, 2, 3, 4] 
     858            sage: a = P._first_smallest_nontrivial(2) 
     859            sage: Integer(a).binary().zfill(32) 
     860            '00000000000000000000000000111100' 
     861            sage: P 
     862            ({0,1,2,3}) 
     863            ({0,1,2,3}) 
     864            ({0,1,2,3}) 
     865            ({0,1,2,3}) 
     866            ({0,1,2,3,4,5}) 
     867            ({0},{1,2,3,4,5}) 
     868            ({0},{1},{2,3,4,5}) 
     869            ({0},{1},{2},{3,4,5}) 
     870            ({0},{1},{2},{3},{4,5}) 
     871            ({0},{1},{2},{3},{4},{5}) 
     872 
     873        """ 
     874        return self.first_smallest_nontrivial(k) 
     875 
     876    cdef unsigned int first_smallest_nontrivial(self, int k): 
     877        cdef unsigned int cell 
     878        cdef int i = 0, j = 0, location = 0, ncols = self.ncols 
     879        cdef int *self_col_lvls = self.col_lvls 
     880        while True: 
     881            if self_col_lvls[i] <= k: 
     882                if i != j and ncols > i - j + 1: 
     883                    ncols = i - j + 1 
     884                    location = j 
     885                j = i + 1 
     886            if self_col_lvls[i] == -1: break 
     887            i += 1 
     888        # location now points to the beginning of the first, smallest, 
     889        # nontrivial cell 
     890        j = location 
     891        while True: 
     892            if self_col_lvls[j] <= k: break 
     893            j += 1 
     894        # j now points to the last element of the cell 
     895        i = self.radix - j - 1                 # the cell is represented in binary, reading from the right: 
     896        cell = (~0 << location) ^ (~0 << j+1)  # <-------            self.radix               -----> 
     897        return cell                            # [0]*(radix-j-1) + [1]*(j-location+1) + [0]*location 
     898 
     899    def _dangerous_dont_use_set_ents_lvls(self, col_ents, col_lvls, wd_ents, wd_lvls): 
     900        """ 
     901        For debugging only. 
     902 
     903        EXAMPLE: 
     904            sage: import sage.coding.binary_code 
     905            sage: from sage.coding.binary_code import * 
     906            sage: P = PartitionStack(4, 6) 
     907            sage: P 
     908            ({0,1,2,3}) 
     909            ({0,1,2,3}) 
     910            ({0,1,2,3}) 
     911            ({0,1,2,3}) 
     912            ({0,1,2,3,4,5}) 
     913            ({0,1,2,3,4,5}) 
     914            ({0,1,2,3,4,5}) 
     915            ({0,1,2,3,4,5}) 
     916            ({0,1,2,3,4,5}) 
     917            ({0,1,2,3,4,5}) 
     918            sage: P._dangerous_dont_use_set_ents_lvls([99]*6, [0,3,2,3,5,-1], [4,3,5,6], [3,2,1,-1]) 
     919            sage: P 
     920            ({4,3,5,6}) 
     921            ({4,3,5},{6}) 
     922            ({4,3},{5},{6}) 
     923            ({4},{3},{5},{6}) 
     924            ({99},{99,99,99,99,99}) 
     925            ({99},{99,99,99,99,99}) 
     926            ({99},{99,99},{99,99,99}) 
     927            ({99},{99},{99},{99},{99,99}) 
     928            ({99},{99},{99},{99},{99,99}) 
     929            ({99},{99},{99},{99},{99},{99}) 
     930 
     931        """ 
     932        cdef unsigned int i 
     933        for i from 0 <= i < len(col_ents): 
     934            self.col_ents[i] = col_ents[i] 
     935        for i from 0 <= i < len(col_lvls): 
     936            self.col_lvls[i] = col_lvls[i] 
     937        for i from 0 <= i < len(wd_ents): 
     938            self.wd_ents[i] = wd_ents[i] 
     939        for i from 0 <= i < len(wd_lvls): 
     940            self.wd_lvls[i] = wd_lvls[i] 
     941 
     942    def _col_percolate(self, start, end): 
     943        """ 
     944        Do one round of bubble sort on ents. 
     945 
     946        EXAMPLE: 
     947            sage: import sage.coding.binary_code 
     948            sage: from sage.coding.binary_code import * 
     949            sage: P = PartitionStack(4, 6) 
     950            sage: P._dangerous_dont_use_set_ents_lvls(range(5,-1,-1), [1,2,2,3,3,-1], range(3,-1,-1), [1,1,2,-1]) 
     951            sage: P 
     952            ({3,2,1,0}) 
     953            ({3},{2},{1,0}) 
     954            ({3},{2},{1},{0}) 
     955            ({3},{2},{1},{0}) 
     956            ({5,4,3,2,1,0}) 
     957            ({5},{4,3,2,1,0}) 
     958            ({5},{4},{3},{2,1,0}) 
     959            ({5},{4},{3},{2},{1},{0}) 
     960            ({5},{4},{3},{2},{1},{0}) 
     961            ({5},{4},{3},{2},{1},{0}) 
     962            sage: P._wd_percolate(0,3) 
     963            sage: P._col_percolate(0,5) 
     964            sage: P 
     965            ({0,3,2,1}) 
     966            ({0},{3},{2,1}) 
     967            ({0},{3},{2},{1}) 
     968            ({0},{3},{2},{1}) 
     969            ({0,5,4,3,2,1}) 
     970            ({0},{5,4,3,2,1}) 
     971            ({0},{5},{4},{3,2,1}) 
     972            ({0},{5},{4},{3},{2},{1}) 
     973            ({0},{5},{4},{3},{2},{1}) 
     974            ({0},{5},{4},{3},{2},{1}) 
     975 
     976        """ 
     977        self.col_percolate(start, end) 
    475978 
    476979    cdef void col_percolate(self, int start, int end): 
    477980        cdef int i, temp 
     981        cdef int *self_col_ents = self.col_ents 
    478982        for i from end >= i > start: 
    479             if self.col_ents[i] < self.col_ents[i-1]: 
     983            if self_col_ents[i] < self_col_ents[i-1]: 
    480984                temp = self.col_ents[i] 
    481                 self.col_ents[i] = self.col_ents[i-1] 
    482                 self.col_ents[i-1] = temp 
    483  
    484     cdef void wd_percolate(self, int start, int end): 
    485         cdef int i, temp 
     985                self_col_ents[i] = self_col_ents[i-1] 
     986                self_col_ents[i-1] = temp 
     987 
     988    def _wd_percolate(self, start, end): 
     989        """ 
     990        Do one round of bubble sort on ents. 
     991         
     992        EXAMPLE: 
     993            sage: import sage.coding.binary_code 
     994            sage: from sage.coding.binary_code import * 
     995            sage: P = PartitionStack(4, 6) 
     996            sage: P._dangerous_dont_use_set_ents_lvls(range(5,-1,-1), [1,2,2,3,3,-1], range(3,-1,-1), [1,1,2,-1]) 
     997            sage: P 
     998            ({3,2,1,0}) 
     999            ({3},{2},{1,0}) 
     1000            ({3},{2},{1},{0}) 
     1001            ({3},{2},{1},{0}) 
     1002            ({5,4,3,2,1,0}) 
     1003            ({5},{4,3,2,1,0}) 
     1004            ({5},{4},{3},{2,1,0}) 
     1005            ({5},{4},{3},{2},{1},{0}) 
     1006            ({5},{4},{3},{2},{1},{0}) 
     1007            ({5},{4},{3},{2},{1},{0}) 
     1008            sage: P._wd_percolate(0,3) 
     1009            sage: P._col_percolate(0,5) 
     1010            sage: P 
     1011            ({0,3,2,1}) 
     1012            ({0},{3},{2,1}) 
     1013            ({0},{3},{2},{1}) 
     1014            ({0},{3},{2},{1}) 
     1015            ({0,5,4,3,2,1}) 
     1016            ({0},{5,4,3,2,1}) 
     1017            ({0},{5},{4},{3,2,1}) 
     1018            ({0},{5},{4},{3},{2},{1}) 
     1019            ({0},{5},{4},{3},{2},{1}) 
     1020            ({0},{5},{4},{3},{2},{1}) 
     1021 
     1022        """ 
     1023        self.wd_percolate(start, end) 
     1024 
     1025    cdef void wd_percolate(self, unsigned int start, unsigned int end): 
     1026        cdef unsigned int i, temp 
     1027        cdef unsigned int *self_wd_ents = self.wd_ents 
    4861028        for i from end >= i > start: 
    487             if self.wd_ents[i] < self.wd_ents[i-1]: 
     1029            if self_wd_ents[i] < self_wd_ents[i-1]: 
    4881030                temp = self.wd_ents[i] 
    489                 self.wd_ents[i] = self.wd_ents[i-1] 
    490                 self.wd_ents[i-1] = temp 
    491  
    492     cdef int split_vertex(self, int col_not_wd, int v, int k): 
     1031                self_wd_ents[i] = self_wd_ents[i-1] 
     1032                self_wd_ents[i-1] = temp 
     1033 
     1034    def _split_column(self, int v, int k): 
     1035        """ 
     1036        Split column v out, placing it before the rest of the cell it was in. 
     1037 
     1038        EXAMPLE: 
     1039            sage: import sage.coding.binary_code 
     1040            sage: from sage.coding.binary_code import * 
     1041            sage: P = PartitionStack(4, 6) 
     1042            sage: [P._split_column(i,i+1) for i in range(5)] 
     1043            [0, 1, 2, 3, 4] 
     1044            sage: P 
     1045            ({0,1,2,3}) 
     1046            ({0,1,2,3}) 
     1047            ({0,1,2,3}) 
     1048            ({0,1,2,3}) 
     1049            ({0,1,2,3,4,5}) 
     1050            ({0},{1,2,3,4,5}) 
     1051            ({0},{1},{2,3,4,5}) 
     1052            ({0},{1},{2},{3,4,5}) 
     1053            ({0},{1},{2},{3},{4,5}) 
     1054            ({0},{1},{2},{3},{4},{5}) 
     1055            sage: P = PartitionStack(4, 6) 
     1056            sage: P._split_column(0,1) 
     1057            0 
     1058            sage: P._split_column(2,2) 
     1059            1 
     1060            sage: P 
     1061            ({0,1,2,3}) 
     1062            ({0,1,2,3}) 
     1063            ({0,1,2,3}) 
     1064            ({0,1,2,3}) 
     1065            ({0,2,1,3,4,5}) 
     1066            ({0},{2,1,3,4,5}) 
     1067            ({0},{2},{1,3,4,5}) 
     1068            ({0},{2},{1,3,4,5}) 
     1069            ({0},{2},{1,3,4,5}) 
     1070            ({0},{2},{1,3,4,5}) 
     1071 
     1072        """ 
     1073        return self.split_column(v, k) 
     1074 
     1075    cdef int split_column(self, int v, int k): 
    4931076        cdef int i = 0, j 
    494         if col_not_wd: 
    495             while self.col_ents[i] != v: i += 1 
    496             j = i 
    497             while self.col_lvls[i] > k: i += 1 
    498             if j == 0 or self.col_lvls[j-1] <= k: 
    499                 self.col_percolate(j+1, i) 
    500             else: 
    501                 while j != 0 and self.col_lvls[j-1] > k: 
    502                     self.col_ents[j] = self.col_ents[j-1] 
    503                     j -= 1 
    504                 self.col_ents[j] = v 
    505             self.col_lvls[j] = k 
    506             return j 
     1077        cdef int *self_col_ents = self.col_ents 
     1078        cdef int *self_col_lvls = self.col_lvls 
     1079        while self_col_ents[i] != v: i += 1 
     1080        j = i 
     1081        while self_col_lvls[i] > k: i += 1 
     1082        if j == 0 or self_col_lvls[j-1] <= k: 
     1083            self.col_percolate(j+1, i) 
    5071084        else: 
    508             while self.wd_ents[i] != v: i += 1 
    509             j = i 
    510             while self.wd_lvls[i] > k: i += 1 
    511             if j == 0 or self.wd_lvls[j-1] <= k: 
    512                 self.wd_percolate(j+1, i) 
    513             else: 
    514                 while j != 0 and self.wd_lvls[j-1] > k: 
    515                     self.wd_ents[j] = self.wd_ents[j-1] 
    516                     j -= 1 
    517                 self.wd_ents[j] = v 
    518             self.wd_lvls[j] = k 
    519             return j 
    520  
    521     cdef int col_degree(self, BinaryCode CG, int col, int wd_ptr, int k): 
    522         cdef int i = 0 
     1085            while j != 0 and self_col_lvls[j-1] > k: 
     1086                self_col_ents[j] = self_col_ents[j-1] 
     1087                j -= 1 
     1088            self_col_ents[j] = v 
     1089        self_col_lvls[j] = k 
     1090        return j 
     1091 
     1092    def _col_degree(self, C, col, wd_ptr, k): 
     1093        """ 
     1094        Returns the number of words in the cell specified by wd_ptr that have a 
     1095        1 in the col-th column. 
     1096 
     1097        EXAMPLE: 
     1098            sage: import sage.coding.binary_code 
     1099            sage: from sage.coding.binary_code import * 
     1100            sage: P = PartitionStack(4, 6) 
     1101            sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]) 
     1102            sage: B = BinaryCode(M) 
     1103            sage: B 
     1104            Binary [6,2] linear code, generator matrix 
     1105            [111100] 
     1106            [001111] 
     1107            sage: [P._split_column(i,i+1) for i in range(5)] 
     1108            [0, 1, 2, 3, 4] 
     1109            sage: P 
     1110            ({0,1,2,3}) 
     1111            ({0,1,2,3}) 
     1112            ({0,1,2,3}) 
     1113            ({0,1,2,3}) 
     1114            ({0,1,2,3,4,5}) 
     1115            ({0},{1,2,3,4,5}) 
     1116            ({0},{1},{2,3,4,5}) 
     1117            ({0},{1},{2},{3,4,5}) 
     1118            ({0},{1},{2},{3},{4,5}) 
     1119            ({0},{1},{2},{3},{4},{5}) 
     1120            sage: P._col_degree(B, 2, 0, 2) 
     1121            2L 
     1122 
     1123        """ 
     1124        return self.col_degree(C, col, wd_ptr, k) 
     1125 
     1126    cdef unsigned int col_degree(self, BinaryCode CG, int col, unsigned int wd_ptr, int k): 
     1127        cdef unsigned int i = 0 
     1128        cdef int *self_wd_lvls = self.wd_lvls 
    5231129        col = self.col_ents[col] 
    5241130        while True: 
    5251131            if CG.is_one(wd_ptr, col): i += 1 
    526             if self.wd_lvls[wd_ptr] > k: wd_ptr += 1 
     1132            if self_wd_lvls[wd_ptr] > k: wd_ptr += 1 
    5271133            else: break 
    5281134        return i 
    5291135 
    530     cdef int wd_degree(self, BinaryCode CG, int wd, int col_ptr, int k): 
     1136    def _wd_degree(self, C, wd, col_ptr, k): 
     1137        """ 
     1138        Returns the number of columns in the cell specified by col_ptr that are 
     1139        1 in wd. 
     1140 
     1141        EXAMPLE: 
     1142            sage: import sage.coding.binary_code 
     1143            sage: from sage.coding.binary_code import * 
     1144            sage: P = PartitionStack(4, 6) 
     1145            sage: M = Matrix(GF(2), [[1,1,1,1,0,0],[0,0,1,1,1,1]]) 
     1146            sage: B = BinaryCode(M) 
     1147            sage: B 
     1148            Binary [6,2] linear code, generator matrix 
     1149            [111100] 
     1150            [001111] 
     1151            sage: [P._split_column(i,i+1) for i in range(5)] 
     1152            [0, 1, 2, 3, 4] 
     1153            sage: P 
     1154            ({0,1,2,3}) 
     1155            ({0,1,2,3}) 
     1156            ({0,1,2,3}) 
     1157            ({0,1,2,3}) 
     1158            ({0,1,2,3,4,5}) 
     1159            ({0},{1,2,3,4,5}) 
     1160            ({0},{1},{2,3,4,5}) 
     1161            ({0},{1},{2},{3,4,5}) 
     1162            ({0},{1},{2},{3},{4,5}) 
     1163            ({0},{1},{2},{3},{4},{5}) 
     1164            sage: P._wd_degree(B, 1, 1, 1) 
     1165            3 
     1166 
     1167        """ 
     1168        return self.wd_degree(C, wd, col_ptr, k) 
     1169 
     1170    cdef int wd_degree(self, BinaryCode CG, unsigned int wd, int col_ptr, int k): 
    5311171        cdef int i = 0 
    5321172        wd = self.wd_ents[wd] 
     1173        cdef int *self_col_lvls = self.col_lvls 
    5331174        while True: 
    5341175            if CG.is_one(wd, col_ptr): i += 1 
    535             if self.col_lvls[col_ptr] > k: col_ptr += 1 
     1176            if self_col_lvls[col_ptr] > k: col_ptr += 1 
    5361177            else: break 
    5371178        return i 
    5381179 
    539     cdef int sort_cols(self, int start, int *degrees, int k): 
    540         cdef int i, j, max, max_location 
    541         cdef int *counts = degrees + self.ncols, *output = degrees + 2*self.ncols 
    542  
    543         for i from 0 <= i < self.ncols: 
    544             counts[i] = 0 
     1180    def _sort_cols(self, start, degrees, k): 
     1181        """ 
     1182        Essentially a counting sort, but on only one cell of the partition. 
     1183 
     1184        INPUT: 
     1185            start -- location of the beginning of the cell 
     1186            k -- at what level of refinement the partition of interest lies 
     1187            degrees -- the counts to sort by 
     1188 
     1189        EXAMPLE: 
     1190            sage: import sage.coding.binary_code 
     1191            sage: from sage.coding.binary_code import * 
     1192            sage: P = PartitionStack(4, 6) 
     1193            sage: [P._split_column(i,i+1) for i in range(5)] 
     1194            [0, 1, 2, 3, 4] 
     1195            sage: P._sort_cols(1, [0,2,2,1,1], 1) 
     1196            2 
     1197            sage: P 
     1198            ({0,1,2,3}) 
     1199            ({0,1,2,3}) 
     1200            ({0,1,2,3}) 
     1201            ({0,1,2,3}) 
     1202            ({0,1,4,5,2,3}) 
     1203            ({0},{1},{4,5},{2,3}) 
     1204            ({0},{1},{4,5},{2,3}) 
     1205            ({0},{1},{4},{5},{2,3}) 
     1206            ({0},{1},{4},{5},{2,3}) 
     1207            ({0},{1},{4},{5},{2},{3}) 
     1208 
     1209        """ 
     1210        cdef int i 
     1211        for i from 0 <= i < len(degrees): 
     1212            self.col_degs[i] = degrees[i] 
     1213        return self.sort_cols(start, k) 
     1214 
     1215    cdef int sort_cols(self, int start, int k): 
     1216        cdef int i, j, max, max_location, self_ncols = self.ncols 
     1217        cdef int *self_col_counts = self.col_counts 
     1218        cdef int *self_col_lvls = self.col_lvls 
     1219        cdef unsigned int *self_col_degs = self.col_degs 
     1220        cdef int *self_col_ents = self.col_ents 
     1221        cdef int *self_col_output = self.col_output 
     1222        for i from 0 <= i < self_ncols: 
     1223            self_col_counts[i] = 0 
    5451224        i = 0 
    546         while self.col_lvls[i+start] > k: 
    547             counts[degrees[i]] += 1 
     1225        while self_col_lvls[i+start] > k: 
     1226            self_col_counts[self_col_degs[i]] += 1 
    5481227            i += 1 
    549         counts[degrees[i]] += 1 
     1228        self_col_counts[self_col_degs[i]] += 1 
    5501229 
    5511230        # i+start is the right endpoint of the cell now 
    552         max = counts[0] 
     1231        max = self_col_counts[0] 
    5531232        max_location = 0 
    554         for j from 0 < j < self.ncols: 
    555             if counts[j] > max: 
    556                 max = counts[j] 
     1233        for j from 0 < j < self_ncols: 
     1234            if self_col_counts[j] > max: 
     1235                max = self_col_counts[j] 
    5571236                max_location = j 
    558             counts[j] += counts[j-1] 
     1237            self_col_counts[j] += self_col_counts[j-1] 
    5591238 
    5601239        for j from i >= j >= 0: 
    561             counts[degrees[j]] -= 1 
    562             output[counts[degrees[j]]] = self.col_ents[start+j] 
    563  
    564         max_location = counts[max_location] + start 
     1240            self_col_counts[self_col_degs[j]] -= 1 
     1241            self_col_output[self_col_counts[self_col_degs[j]]] = self_col_ents[start+j] 
     1242 
     1243        max_location = self_col_counts[max_location] + start 
    5651244 
    5661245        for j from 0 <= j <= i: 
    567             self.col_ents[start+j] = output[j] 
     1246            self_col_ents[start+j] = self_col_output[j] 
    5681247 
    5691248        j = 1 
    570         while j < self.ncols and counts[j] <= i: 
    571             if counts[j] > 0: 
    572                 self.col_lvls[start + counts[j] - 1] = k 
    573             self.col_percolate(start + counts[j-1], start + counts[j] - 1) 
     1249        while j < self_ncols and self_col_counts[j] <= i: 
     1250            if self_col_counts[j] > 0: 
     1251                self_col_lvls[start + self_col_counts[j] - 1] = k 
     1252            self.col_percolate(start + self_col_counts[j-1], start + self_col_counts[j] - 1) 
    5741253            j += 1 
    5751254 
    5761255        return max_location 
    5771256 
    578     cdef int sort_wds(self, int start, int *degrees, int k): 
    579         cdef int i, j, max, max_location 
    580         cdef int *counts = degrees + self.nwords, *output = degrees + 2*self.nwords 
    581  
    582         for i from 0 <= i < self.nwords: 
    583             counts[i] = 0 
     1257    def _sort_wds(self, start, degrees, k): 
     1258        """ 
     1259        Essentially a counting sort, but on only one cell of the partition. 
     1260 
     1261        INPUT: 
     1262            start -- location of the beginning of the cell 
     1263            k -- at what level of refinement the partition of interest lies 
     1264            degrees -- the counts to sort by 
     1265 
     1266        EXAMPLE: 
     1267            sage: import sage.coding.binary_code 
     1268            sage: from sage.coding.binary_code import * 
     1269            sage: P = PartitionStack(8, 6) 
     1270            sage: P._sort_wds(0, [0,0,3,3,3,3,2,2], 1) 
     1271            4L 
     1272            sage: P 
     1273            ({0,1,6,7,2,3,4,5}) 
     1274            ({0,1},{6,7},{2,3,4,5}) 
     1275            ({0,1},{6,7},{2,3,4,5}) 
     1276            ({0,1},{6,7},{2,3,4,5}) 
     1277            ({0,1},{6,7},{2,3,4,5}) 
     1278            ({0,1},{6,7},{2,3,4,5}) 
     1279            ({0,1},{6,7},{2,3,4,5}) 
     1280            ({0,1},{6,7},{2,3,4,5}) 
     1281            ({0,1,2,3,4,5}) 
     1282            ({0,1,2,3,4,5}) 
     1283            ({0,1,2,3,4,5}) 
     1284            ({0,1,2,3,4,5}) 
     1285            ({0,1,2,3,4,5}) 
     1286            ({0,1,2,3,4,5}) 
     1287 
     1288        """ 
     1289        cdef int i 
     1290        for i from 0 <= i < len(degrees): 
     1291            self.wd_degs[i] = degrees[i] 
     1292        return self.sort_wds(start, k) 
     1293 
     1294    cdef unsigned int sort_wds(self, unsigned int start, int k): 
     1295        cdef unsigned int i, j, max, max_location, self_nwords = self.nwords 
     1296        cdef unsigned int *self_wd_counts = self.wd_counts 
     1297        cdef int *self_wd_lvls = self.wd_lvls 
     1298        cdef int *self_wd_degs = self.wd_degs 
     1299        cdef unsigned int *self_wd_ents = self.wd_ents 
     1300        cdef unsigned int *self_wd_output = self.wd_output 
     1301 
     1302        for i from 0 <= i < self_nwords: 
     1303            self_wd_counts[i] = 0 
    5841304        i = 0 
    585         while self.wd_lvls[i+start] > k: 
    586             counts[degrees[i]] += 1 
     1305        while self_wd_lvls[i+start] > k: 
     1306            self_wd_counts[self_wd_degs[i]] += 1 
    5871307            i += 1 
    588         counts[degrees[i]] += 1 
     1308        self_wd_counts[self_wd_degs[i]] += 1 
    5891309 
    5901310        # i+start is the right endpoint of the cell now 
    591         max = counts[0] 
     1311        max = self_wd_counts[0] 
    5921312        max_location = 0 
    593         for j from 0 < j < self.nwords: 
    594             if counts[j] > max: 
    595                 max = counts[j] 
     1313        for j from 0 < j < self_nwords: 
     1314            if self_wd_counts[j] > max: 
     1315                max = self_wd_counts[j] 
    5961316                max_location = j 
    597             counts[j] += counts[j-1] 
     1317            self_wd_counts[j] += self_wd_counts[j-1] 
    5981318 
    5991319        for j from i >= j >= 0: 
    600             counts[degrees[j]] -= 1 
    601             output[counts[degrees[j]]] = self.wd_ents[start+j] 
    602  
    603         max_location = counts[max_location] + start 
     1320            if j > i: break # cython bug with unsigned ints... 
     1321            self_wd_counts[self_wd_degs[j]] -= 1 
     1322            self_wd_output[self_wd_counts[self_wd_degs[j]]] = self_wd_ents[start+j] 
     1323 
     1324        max_location = self_wd_counts[max_location] + start 
    6041325 
    6051326        for j from 0 <= j <= i: 
    606             self.wd_ents[start+j] = output[j] 
     1327            self_wd_ents[start+j] = self_wd_output[j] 
    6071328 
    6081329        j = 1 
    609         while j < self.nwords and counts[j] <= i: 
    610             if counts[j] > 0: 
    611                 self.wd_lvls[start + counts[j] - 1] = k 
    612             self.wd_percolate(start + counts[j-1], start + counts[j] - 1) 
     1330        while j < self_nwords and self_wd_counts[j] <= i: 
     1331            if self_wd_counts[j] > 0: 
     1332                self_wd_lvls[start + self_wd_counts[j] - 1] = k 
     1333            self.wd_percolate(start + self_wd_counts[j-1], start + self_wd_counts[j] - 1) 
    6131334            j += 1 
    6141335 
    6151336        return max_location 
    6161337 
    617     cdef int refine(self, int k, int *col_alpha, int *wd_alpha, BinaryCode CG): 
     1338################################################################################ 
     1339################################################################################ 
     1340################################################################################ 
     1341 
     1342    def _refine(self, k, col_alpha, wd_alpha, Code): 
     1343        """ 
     1344        Refines the partition at level k, using the list of cells alpha, and Code. 
     1345 
     1346        EXAMPLE: 
     1347 
     1348        """ 
     1349        cdef unsigned int i 
     1350        cdef int j 
     1351        cdef unsigned int *_col_a = <unsigned int *> sage_malloc(4*self.ncols*sizeof(unsigned int)) 
     1352        cdef int *_wd_a = <int *> sage_malloc(4*self.nwords*sizeof(int)) 
     1353        if not _col_a or not _wd_a: 
     1354            if _col_a: sage_free(_col_a) 
     1355            if _wd_a: sage_free(_wd_a) 
     1356            raise MemoryError("Memory.") 
     1357    # TODO       for i from  TODO 
     1358        result = self.refine(k, _col_a, _wd_a, Code) 
     1359        sage_free(_col_a) 
     1360        sage_free(_wd_a) 
     1361        return result 
     1362 
     1363    cdef unsigned int refine(self, int k, unsigned int *col_alpha, int *wd_alpha, BinaryCode CG): 
    6181364        cdef int m = 0, j 
    6191365        cdef int i, q, r, s, t 
    620         cdef int invariant 
    621         cdef int *col_degrees = col_alpha + self.ncols 
     1366        cdef unsigned int t_w 
     1367        cdef unsigned int invariant 
     1368        cdef unsigned int *col_degrees = col_alpha + self.ncols 
    6221369        cdef int *wd_degrees = wd_alpha + self.nwords 
    6231370        while not self.is_discrete(k) and col_alpha[m] != -1: 
    … …  
    6341381                if s: 
    6351382                    invariant += 8 
    636                     t = self.sort_wds(j, wd_degrees, k) 
    637                     invariant += t + wd_degrees[i-j-1] 
     1383#FIX                 t_w = self.sort_wds(j, wd_degrees, k) 
     1384                    invariant += t_w + wd_degrees[i-j-1] 
    6381385                    q = m 
    6391386                    while col_alpha[q] != -1: 
    640                         if col_alpha[q] == j: col_alpha[q] = t 
     1387                        if col_alpha[q] == j: col_alpha[q] = t_w 
    6411388                        q += 1 
    6421389                    r = j 
    6431390                    while True: 
    6441391                        if r == 0 or self.wd_lvls[r-1] == k: 
    645                             if r != t: 
     1392                            if r != t_w: 
    6461393                                col_alpha[q] = r 
    6471394                                q += 1 
    … …  
    6691416                if s: 
    6701417                    invariant += 8 
    671                     t = self.sort_cols(j, col_degrees, k) 
     1418#FIX                 t = self.sort_cols(j, col_degrees, k) 
    6721419                    invariant += t + col_degrees[i-j-1] 
    6731420                    q = m 
    … …  
    7001447 
    7011448    cdef int cmp(self, PartitionStack other, BinaryCode CG): 
    702         # if CG(self) > G(other): return 1 
    703         # if CG(self) < G(other): return -1 
     1449        # if CG(self) > CG(other): return 1 
     1450        # if CG(self) < CG(other): return -1 
    7041451        # else: return 0 
    7051452        cdef int i, j, k, l 
    … …  
    7081455                k = CG.is_one(self.wd_ents[i], self.col_ents[j]) 
    7091456                l = CG.is_one(other.wd_ents[i], other.col_ents[j]) 
    710                 if k ^ l: 
     1457                if k != l: 
    7111458                    return k - l 
    7121459        return 0 
     1460 
     1461cdef class BinaryCodeClassifier: 
     1462 
     1463    def __new__(self): 
     1464        cdef int *self_ham_wts 
     1465        self.ham_wts = <int *> sage_malloc( 65536 * sizeof(int) ) 
     1466        if not self.ham_wts: 
     1467            sage_free(self.ham_wts) 
     1468        self_ham_wts = self.ham_wts 
     1469        self_ham_wts[0]  = 0; self_ham_wts[1]  = 1; self_ham_wts[2]  = 1; self_ham_wts[3]  = 2 
     1470        self_ham_wts[4]  = 1; self_ham_wts[5]  = 2; self_ham_wts[6]  = 2; self_ham_wts[7]  = 3 
     1471        self_ham_wts[8]  = 1; self_ham_wts[9]  = 2; self_ham_wts[10] = 2; self_ham_wts[11] = 3 
     1472        self_ham_wts[12] = 2; self_ham_wts[13] = 3; self_ham_wts[14] = 3; self_ham_wts[15] = 4 
     1473        for i from 16 <= i < 256: 
     1474            self_ham_wts[i] = self_ham_wts[i & 15] + self_ham_wts[(i>>4) & 15] 
     1475        for i from 256 <= i < 65536: 
     1476            self_ham_wts[i] = self_ham_wts[i & 255] + self_ham_wts[(i>>8) & 255] 
     1477        # This may seem like overkill, but the worst case for storing the words 
     1478        # themselves is 65536- in this case, we are increasing memory usage by a 
     1479        # factor of 2. 
     1480 
     1481    def __dealloc__(self): 
     1482        sage_free(self.ham_wts) 
     1483 
     1484 
     1485 
     1486 
     1487 
     1488 
     1489 
     1490 
     1491 
     1492 
     1493 
     1494 
     1495 
     1496 
     1497 
     1498 
     1499 
     1500 
    7131501 
    7141502def classify(BinaryCode C, lab=True, verbosity=0): 

  • sage/graphs/graph_isom.pyx

    r7757 r7760  
    577577                    s = m 
    578578                    while alpha[s] != -1: 
    579                         if alpha[s] == j: alpha[s] = t 
     579                        if alpha[s] == j: alpha[s] = t # TODO this will only happen once, so should break 
    580580                        s += 1 
    581581                    r = j 
    … …  
    621621                    s = m 
    622622                    while alpha[s] != -1: 
    623                         if alpha[s] == j: alpha[s] = t 
     623                        if alpha[s] == j: alpha[s] = t # this will only happen once, so should break 
    624624                        s += 1 
    625625                    r = j 
Note: See TracChangeset for help on using the changeset viewer.