Changeset 7758:aff79c4f15e7


Ignore:
Timestamp:
11/08/07 10:40:20 (6 years ago)
Author:
Robert L Miller <rlm@…>
Branch:
default
Message:

2nd half

Location:
sage/coding
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • sage/coding/binary_code.pxd

    r7756 r7758  
    11 
    2 cdef class BinaryCodeGraph: 
    3      
    4     cdef int *columns 
     2include '../ext/cdefs.pxi' 
     3 
     4cdef class BinaryCode: 
     5    cdef unsigned int *columns 
    56    cdef int *ham_wts 
    67    cdef int ncols 
    78    cdef int nrows 
    89    cdef int radix 
    9     cdef int ptn_mask 
    10     cdef int nwords 
    11     cdef int has_edge_bip(self, int, int) 
    12     cdef int has_edge(self, int, int) 
     10    cdef unsigned int nwords 
     11    cdef int is_one(self, unsigned int, int) 
     12    cdef int is_automorphism(self, int *, unsigned int *) 
     13 
     14cdef class OrbitPartition: 
     15    cdef unsigned int *wd_parent 
     16    cdef unsigned int *wd_rank 
     17    cdef unsigned int *wd_min_cell_rep 
     18    cdef unsigned int *wd_size 
     19    cdef int *col_parent 
     20    cdef int *col_rank 
     21    cdef int *col_min_cell_rep 
     22    cdef int *col_size 
     23    cdef unsigned int wd_find(self, unsigned int) 
     24    cdef void wd_union(self, unsigned int, unsigned int) 
     25    cdef int col_find(self, int) 
     26    cdef void col_union(self, int, int) 
     27    cdef int merge_perm(self, int *, unsigned int *, int, int) 
    1328 
    1429cdef class PartitionStack: 
    1530    cdef int *wd_ents 
     31    cdef int *wd_lvls 
    1632    cdef int *col_ents 
    17     cdef int *wd_lvls 
    1833    cdef int *col_lvls 
    1934    cdef int nwords 
    … …  
    2742    cdef void wd_percolate(self, int start, int end) 
    2843    cdef int split_vertex(self, int, int, int) 
    29     cdef int col_degree(self, BinaryCodeGraph, int, int, int) 
    30     cdef int wd_degree(self, BinaryCodeGraph, int, int, int) 
     44    cdef int col_degree(self, BinaryCode, int, int, int) 
     45    cdef int wd_degree(self, BinaryCode, int, int, int) 
    3146    cdef int sort_cols(self, int, int *, int) 
    3247    cdef int sort_wds(self, int, int *, int) 
    33     cdef int refine(self, int, int *, int *, BinaryCodeGraph) 
     48    cdef int refine(self, int, int *, int *, BinaryCode) 
    3449    cdef void get_permutation(self, PartitionStack, int *, int *) 
    35     cdef int cmp(self, PartitionStack, BinaryCodeGraph) 
     50    cdef int cmp(self, PartitionStack, BinaryCode) 

  • sage/coding/binary_code.pyx

    r7756 r7758  
     1""" 
     2Fast binary code routines. 
     3 
     4Some computations with linear binary codes. Fix a basis for $GF(2)^n$. 
     5A linear binary code is a linear subspace of $GF(2)^n$, together with 
     6this choice of basis. A permutation $g \in S_n$ of the fixed basis 
     7gives rise to a permutation of the vectors, or words, in $GF(2)^n$, 
     8sending $(w_i)$ to $(w_{g(i)})$. The permutation automorphism group of 
     9the code $C$ is the set of permutations of the basis that bijectively 
     10map $C$ to itself. Note that if $g$ is such a permutation, then  
     11$$g(a_i) + g(b_i) = (a_{g(i)} + b_{g(i)}) = g((a_i) + (b_i)).$$ 
     12Over other fields, it is also required that the map be linear, which 
     13as per above boils down to scalar multiplication. However, over 
     14$GF(2),$ the only scalars are 0 and 1, so the linearity condition has 
     15trivial effect. 
     16 
     17AUTHOR: 
     18    Robert L Miller (Oct-Nov 2007) 
     19        * Compiled code datastructure 
     20        * union-find based orbit partition 
     21        * optimized partition stack class 
     22 
     23""" 
     24 
     25#***************************************************************************** 
     26#         Copyright (C) 2007 Robert L. Miller <rlmillster@gmail.com> 
     27# 
     28# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL) 
     29#                         http://www.gnu.org/licenses/ 
     30#***************************************************************************** 
    131 
    232include '../ext/cdefs.pxi' 
    … …  
    434include '../ext/stdsage.pxi' 
    535from sage.structure.element import is_Matrix 
    6 #from sage.graphs.graph_isom cimport OrbitPartition,\ 
    7 #    _orbit_partition_from_list_perm 
    836from sage.misc.misc import cputime 
    937from math import log, ceil 
    1038from sage.rings.integer import Integer 
    1139 
    12 cdef class BinaryCodeGraph: 
    13      
     40cdef class BinaryCode: 
     41    """ 
     42    Minimal, but optimized, binary code object. 
     43 
     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 
     62    """ 
    1463    def __new__(self, arg1): 
    15         cdef unsigned int i, k, j, z         
    16         self.radix = 8*sizeof(int) 
     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 
     83        self.radix = 8*sizeof(unsigned int) 
    1784        self.ncols = arg1.ncols() 
    1885        self.nrows = arg1.nrows() 
    1986        self.nwords = 1 << self.nrows 
    20         self.ptn_mask = ~0 << self.nrows 
    21          
    22         if self.nrows >= self.radix or ceil(log(self.ncols,2)) >= self.radix: 
     87 
     88        if self.nrows > self.radix or ceil(log(self.ncols,2)) > self.radix: 
    2389            raise NotImplementedError("Columns and rows are stored as ints. This code is too big.") 
    24         self.columns = <int *> sage_malloc( self.ncols * sizeof(int) ) 
     90        self.columns = <unsigned int *> sage_malloc( self.ncols * sizeof(unsigned int) ) 
    2591        if not self.columns: 
    2692            raise MemoryError("Memory.") 
    27          
     93 
    2894        cols = arg1.columns() 
    2995        for i from 0 <= i < self.ncols: 
    … …  
    3298                k += (1 << j) 
    3399            self.columns[i] = k 
    34          
     100 
    35101        self.ham_wts = <int *> sage_malloc( 16 * sizeof(int) ) 
    36102        if not self.ham_wts: 
    … …  
    41107        self.ham_wts[8]  = 1; self.ham_wts[9]  = 2; self.ham_wts[10] = 2; self.ham_wts[11] = 3 
    42108        self.ham_wts[12] = 2; self.ham_wts[13] = 3; self.ham_wts[14] = 3; self.ham_wts[15] = 4 
    43                  
     109 
    44110    def __dealloc__(self): 
    45111        sage_free(self.columns) 
    46112        sage_free(self.ham_wts) 
    47113 
    48     cdef int has_edge_bip(self, int word, int column): 
    49         cdef int i, j, k 
     114    def __repr__(self): 
     115        """ 
     116        String representation of self. 
     117 
     118        EXAMPLE: 
     119            sage: import sage.coding.binary_code 
     120            sage: from sage.coding.binary_code import * 
     121            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]]) 
     122            sage: B = BinaryCode(M) 
     123            sage: B 
     124            Binary [8,4] linear code, generator matrix 
     125            [11110000] 
     126            [00111100] 
     127            [00001111] 
     128            [10101010] 
     129 
     130        """ 
     131        cdef int i, j 
     132        s = 'Binary [%d,%d] linear code, generator matrix\n'%(self.ncols, self.nrows) 
     133        for i from 0 <= i < self.nrows: 
     134            s += '[' 
     135            for j from 0 <= j < self.ncols: 
     136                s += '%d'%self.is_one(1<<i,j) 
     137            s += ']\n' 
     138        return s 
     139 
     140    def _is_one(self, word, col): 
     141        """ 
     142        Returns the col-th letter of word, i.e. 0 or 1. Words are expressed 
     143        as integers, which represent linear combinations of the rows of the 
     144        generator matrix of the code. 
     145 
     146        EXAMPLE: 
     147            sage: import sage.coding.binary_code 
     148            sage: from sage.coding.binary_code import * 
     149            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]]) 
     150            sage: B = BinaryCode(M) 
     151            sage: B 
     152            Binary [8,4] linear code, generator matrix 
     153            [11110000] 
     154            [00111100] 
     155            [00001111] 
     156            [10101010] 
     157            sage: B._is_one(7, 4) 
     158            0 
     159            sage: B._is_one(8, 4) 
     160            1 
     161            sage: B._is_automorphism([1,0,3,2,4,5,6,7], [1, 2, 4, 9]) 
     162            1 
     163 
     164        """ 
     165        return self.is_one(word, col) 
     166 
     167    cdef int is_one(self, unsigned int word, int column): 
     168        cdef int i, j 
     169        cdef unsigned int k 
    50170        i = 0 
    51171        k = word & self.columns[column] 
    … …  
    55175        return i % 2 
    56176 
    57     cdef int has_edge(self, int a, int b): 
     177    def _is_automorphism(self, col_gamma, basis_gamma): 
     178        """ 
     179        Check whether a given permutation is an automorphism of the code. 
     180 
     181        INPUT: 
     182            col_gamma -- permutation sending i |--> col_gamma[i] acting 
     183                on the columns. 
     184            basis_gamma -- describes where the basis elements are mapped 
     185                under gamma. basis_gamma[i] is where the i-th row is sent, 
     186                as an integer expressing a linear combination of the rows 
     187                of the generator matrix. 
     188 
     189        EXAMPLE: 
     190            sage: import sage.coding.binary_code 
     191            sage: from sage.coding.binary_code import * 
     192            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]]) 
     193            sage: B = BinaryCode(M) 
     194            sage: B 
     195            Binary [8,4] linear code, generator matrix 
     196            [11110000] 
     197            [00111100] 
     198            [00001111] 
     199            [10101010] 
     200            sage: B._is_automorphism([1,0,3,2,4,5,6,7], [1, 2, 4, 9]) 
     201            1 
     202 
     203        """ 
     204        cdef int i 
     205        cdef int *_col_gamma 
     206        cdef unsigned int *_basis_gamma 
     207        _basis_gamma = <unsigned int *> sage_malloc(self.nrows * sizeof(unsigned int)) 
     208        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int)) 
     209        if not (_col_gamma and _basis_gamma): 
     210            if _basis_gamma: sage_free(_basis_gamma) 
     211            if _col_gamma: sage_free(_col_gamma) 
     212            raise MemoryError("Memory.") 
     213        for i from 0 <= i < self.nrows: 
     214            _basis_gamma[i] = basis_gamma[i] 
     215        for i from 0 <= i < self.ncols: 
     216            _col_gamma[i] = col_gamma[i] 
     217        result = self.is_automorphism(_col_gamma, _basis_gamma) 
     218        sage_free(_col_gamma) 
     219        sage_free(_basis_gamma) 
     220        return result 
     221 
     222    cdef int is_automorphism(self, int *col_gamma, unsigned int *basis_gamma): 
    58223        cdef int i, j 
    59         i = self.ptn_mask & a 
    60         j = self.ptn_mask & b 
    61         if (i==0)==(0==j): 
    62             # both are on the same side of the bipartite graph 
    63             return 0 
    64         elif i: 
    65             # first argument is a column 
    66             return self.has_edge_bip(b, a - self.nwords) 
     224        for i from 0 <= i < self.nrows: 
     225            for j from 0 <= j < self.ncols: 
     226                if self.is_one(1 << i, j) != self.is_one(basis_gamma[i], col_gamma[j]): 
     227                    return 0 
     228        return 1 
     229 
     230cdef class OrbitPartition: 
     231    """ 
     232    Structure which keeps track of which vertices are equivalent 
     233    under the part of the automorphism group that has already been 
     234    seen, during search. Essentially a disjoint-set data structure*, 
     235    which also keeps track of the minimum element and size of each 
     236    cell of the partition, and the size of the partition. 
     237 
     238    * http://en.wikipedia.org/wiki/Disjoint-set_data_structure 
     239 
     240    """ 
     241    def __new__(self, nrows, ncols): 
     242        cdef unsigned int nwords, word 
     243        cdef int col 
     244        nwords = (1 << nrows) 
     245        self.wd_parent =       <unsigned int *> sage_malloc( nwords * sizeof(unsigned int) ) 
     246        self.wd_rank =         <unsigned int *> sage_malloc( nwords * sizeof(unsigned int) ) 
     247        self.wd_min_cell_rep = <unsigned int *> sage_malloc( nwords * sizeof(unsigned int) ) 
     248        self.wd_size =         <unsigned int *> sage_malloc( nwords * sizeof(unsigned int) ) 
     249        self.col_parent =       <int *> sage_malloc( ncols * sizeof(int) ) 
     250        self.col_rank =         <int *> sage_malloc( ncols * sizeof(int) ) 
     251        self.col_min_cell_rep = <int *> sage_malloc( ncols * sizeof(int) ) 
     252        self.col_size =         <int *> sage_malloc( ncols * sizeof(int) ) 
     253        if not (self.wd_parent and self.wd_rank and self.wd_min_cell_rep and self.wd_size and self.col_parent and self.col_rank and self.col_min_cell_rep and self.col_size): 
     254            if self.wd_parent: sage_free(self.wd_parent) 
     255            if self.wd_rank: sage_free(self.wd_rank) 
     256            if self.wd_min_cell_rep: sage_free(self.wd_min_cell_rep) 
     257            if self.wd_size: sage_free(self.wd_size) 
     258            if self.col_parent: sage_free(self.col_parent) 
     259            if self.col_rank: sage_free(self.col_rank) 
     260            if self.col_min_cell_rep: sage_free(self.col_min_cell_rep) 
     261            if self.col_size: sage_free(self.col_size) 
     262            raise MemoryError("Memory.") 
     263        for word from 0 <= word < nwords: 
     264            self.wd_parent[word] = word 
     265            self.wd_rank[word] = 0 
     266            self.wd_min_cell_rep[word] = word 
     267            self.wd_size[word] = 1 
     268        for col from 0 <= col < ncols: 
     269            self.col_parent[col] = col 
     270            self.col_rank[col] = 0 
     271            self.col_min_cell_rep[col] = col 
     272            self.col_size[col] = 1 
     273 
     274    def __dealloc__(self): 
     275        sage_free(self.wd_parent) 
     276        sage_free(self.wd_rank) 
     277        sage_free(self.wd_min_cell_rep) 
     278        sage_free(self.wd_size) 
     279        sage_free(self.col_parent) 
     280        sage_free(self.col_rank) 
     281        sage_free(self.col_min_cell_rep) 
     282        sage_free(self.col_size) 
     283 
     284    cdef unsigned int wd_find(self, unsigned int word): 
     285        if self.wd_parent[word] == word: 
     286            return word 
    67287        else: 
    68             # first argument is a vector 
    69             return self.has_edge_bip(a, b - self.nwords) 
    70      
    71 #    cdef int is_automorphism(self, int *wd_gamma, int *col_gamma): 
    72          
    73         #???? 
    74          
    75         #cdef int i, j, k, l 
    76         #for i from 0 <= i < self.nwords: 
    77         #    for j from 0 <= j < self.ncols: 
    78         #        k = CG.has_edge_bip(self.wd_ents[i], self.col_ents[j]) 
    79         #        l = CG.has_edge_bip(other.wd_ents[i], other.col_ents[j]) 
    80         #        if k ^ l: 
    81         #            return k - l 
    82         #return 0 
    83      
     288            self.wd_parent[word] = self.wd_find(self.wd_parent[word]) 
     289            return self.wd_parent[word] 
     290 
     291    cdef void wd_union(self, unsigned int x, unsigned int y): 
     292        cdef unsigned int x_root, y_root 
     293        x_root = self.wd_find(x) 
     294        y_root = self.wd_find(y) 
     295        if self.wd_rank[x_root] > self.wd_rank[y_root]: 
     296            self.wd_parent[y_root] = x_root 
     297            self.wd_min_cell_rep[y_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root]) 
     298            self.wd_size[y_root] += self.wd_size[x_root] 
     299        elif self.wd_rank[x_root] < self.wd_rank[y_root]: 
     300            self.wd_parent[x_root] = y_root 
     301            self.wd_min_cell_rep[x_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root]) 
     302            self.wd_size[x_root] += self.wd_size[y_root] 
     303        elif x_root != y_root: 
     304            self.wd_parent[y_root] = x_root 
     305            self.wd_min_cell_rep[y_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root]) 
     306            self.wd_size[y_root] += self.wd_size[x_root] 
     307            self.wd_rank[x_root] += 1 
     308 
     309    cdef int col_find(self, int col): 
     310        if self.col_parent[col] == col: 
     311            return col 
     312        else: 
     313            self.col_parent[col] = self.col_find(self.col_parent[col]) 
     314            return self.col_parent[col] 
     315 
     316    cdef void col_union(self, int x, int y): 
     317        cdef int x_root, y_root 
     318        x_root = self.col_find(x) 
     319        y_root = self.col_find(y) 
     320        if self.col_rank[x_root] > self.col_rank[y_root]: 
     321            self.col_parent[y_root] = x_root 
     322            self.col_min_cell_rep[y_root] = min(self.col_min_cell_rep[x_root],self.col_min_cell_rep[y_root]) 
     323            self.col_size[y_root] += self.col_size[x_root] 
     324        elif self.col_rank[x_root] < self.col_rank[y_root]: 
     325            self.col_parent[x_root] = y_root 
     326            self.col_min_cell_rep[x_root] = min(self.col_min_cell_rep[x_root],self.col_min_cell_rep[y_root]) 
     327            self.col_size[x_root] += self.col_size[y_root] 
     328        elif x_root != y_root: 
     329            self.col_parent[y_root] = x_root 
     330            self.col_min_cell_rep[y_root] = min(self.col_min_cell_rep[x_root],self.col_min_cell_rep[y_root]) 
     331            self.col_size[y_root] += self.col_size[x_root] 
     332            self.col_rank[x_root] += 1 
     333 
     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        """ 
     341        cdef unsigned int i, gamma_i_root 
     342        cdef int j, gamma_j_root, return_value = 0 
     343        for i from 0 <= i < (1 << nrows): 
     344            if self.wd_parent[i] == i: 
     345                gamma_i_root = self.wd_find(wd_gamma[i]) 
     346                if gamma_i_root != i: 
     347                    return_value = 1 
     348                    self.wd_union(i, gamma_i_root) 
     349        for j from 0 <= j < ncols: 
     350            if self.col_parent[j] == j: 
     351                gamma_j_root = self.col_find(col_gamma[j]) 
     352                if gamma_j_root != j: 
     353                    return_value = 1 
     354                    self.col_union(j, gamma_j_root) 
     355        return return_value 
     356 
    84357cdef class PartitionStack: 
    85      
     358 
    86359    def __new__(self, arg1, arg2=None): 
    87360        cdef int k 
    … …  
    97370            if self.col_ents: sage_free(self.col_ents) 
    98371            if self.col_lvls: sage_free(self.col_lvls) 
     372            raise MemoryError("Memory.") 
    99373        for k from 0 <= k < self.nwords-1: 
    100374            self.wd_ents[k] = k 
    … …  
    107381        self.col_ents[self.ncols-1] = self.ncols-1 
    108382        self.col_lvls[self.ncols-1] = -1 
    109      
     383 
    110384    def __dealloc__(self): 
    111385        sage_free(self.wd_ents) 
    … …  
    136410            s = s[:-2] + ')\n' 
    137411        return s 
    138          
     412 
    139413    cdef int is_discrete(self, int k): 
    140414        cdef int i 
    … …  
    143417                return 0 
    144418        return 1 
    145      
     419 
    146420    cdef int num_cells(self, int k): 
    147421        cdef int i, j = 0 
    … …  
    153427                j += 1 
    154428        return j 
    155      
     429 
    156430    cdef int sat_225(self, int k): 
    157431        cdef int i, n = self.nwords + self.ncols, in_cell = 0 
    … …  
    179453            return 1 
    180454        return 0 
    181          
     455 
    182456    cdef int is_min_cell_rep(self, int col_not_wd, int i, int k): 
    183457        if i == 0: 
    … …  
    187461        else: 
    188462            return self.wd_lvls[i-1] <= k 
    189      
     463 
    190464    cdef int is_fixed(self, int col_not_wd, int i, int k): 
    191465        """ 
    … …  
    196470        else: 
    197471            return self.wd_lvls[i] <= k 
    198      
     472 
    199473    ## TODO: first smallest nontrivial 
    200      
     474 
     475 
    201476    cdef void col_percolate(self, int start, int end): 
    202477        cdef int i, temp 
    … …  
    206481                self.col_ents[i] = self.col_ents[i-1] 
    207482                self.col_ents[i-1] = temp 
    208      
     483 
    209484    cdef void wd_percolate(self, int start, int end): 
    210485        cdef int i, temp 
    … …  
    214489                self.wd_ents[i] = self.wd_ents[i-1] 
    215490                self.wd_ents[i-1] = temp 
    216      
     491 
    217492    cdef int split_vertex(self, int col_not_wd, int v, int k): 
    218493        cdef int i = 0, j 
    … …  
    243518            self.wd_lvls[j] = k 
    244519            return j 
    245      
    246     cdef int col_degree(self, BinaryCodeGraph CG, int col, int wd_ptr, int k): 
     520 
     521    cdef int col_degree(self, BinaryCode CG, int col, int wd_ptr, int k): 
    247522        cdef int i = 0 
    248523        col = self.col_ents[col] 
    249524        while True: 
    250             if CG.has_edge_bip(wd_ptr, col): i += 1 
     525            if CG.is_one(wd_ptr, col): i += 1 
    251526            if self.wd_lvls[wd_ptr] > k: wd_ptr += 1 
    252527            else: break 
    253528        return i 
    254      
    255     cdef int wd_degree(self, BinaryCodeGraph CG, int wd, int col_ptr, int k): 
     529 
     530    cdef int wd_degree(self, BinaryCode CG, int wd, int col_ptr, int k): 
    256531        cdef int i = 0 
    257532        wd = self.wd_ents[wd] 
    258533        while True: 
    259             if CG.has_edge_bip(wd, col_ptr): i += 1 
     534            if CG.is_one(wd, col_ptr): i += 1 
    260535            if self.col_lvls[col_ptr] > k: col_ptr += 1 
    261536            else: break 
    262537        return i 
    263      
     538 
    264539    cdef int sort_cols(self, int start, int *degrees, int k): 
    265540        cdef int i, j, max, max_location 
    266541        cdef int *counts = degrees + self.ncols, *output = degrees + 2*self.ncols 
    267          
     542 
    268543        for i from 0 <= i < self.ncols: 
    269544            counts[i] = 0 
    … …  
    273548            i += 1 
    274549        counts[degrees[i]] += 1 
    275          
     550 
    276551        # i+start is the right endpoint of the cell now 
    277552        max = counts[0] 
    … …  
    282557                max_location = j 
    283558            counts[j] += counts[j-1] 
    284          
     559 
    285560        for j from i >= j >= 0: 
    286561            counts[degrees[j]] -= 1 
    287562            output[counts[degrees[j]]] = self.col_ents[start+j] 
    288          
     563 
    289564        max_location = counts[max_location] + start 
    290          
     565 
    291566        for j from 0 <= j <= i: 
    292567            self.col_ents[start+j] = output[j] 
    293          
     568 
    294569        j = 1 
    295570        while j < self.ncols and counts[j] <= i: 
    … …  
    298573            self.col_percolate(start + counts[j-1], start + counts[j] - 1) 
    299574            j += 1 
    300          
     575 
    301576        return max_location 
    302          
     577 
    303578    cdef int sort_wds(self, int start, int *degrees, int k): 
    304579        cdef int i, j, max, max_location 
    305580        cdef int *counts = degrees + self.nwords, *output = degrees + 2*self.nwords 
    306          
     581 
    307582        for i from 0 <= i < self.nwords: 
    308583            counts[i] = 0 
    … …  
    312587            i += 1 
    313588        counts[degrees[i]] += 1 
    314          
     589 
    315590        # i+start is the right endpoint of the cell now 
    316591        max = counts[0] 
    … …  
    321596                max_location = j 
    322597            counts[j] += counts[j-1] 
    323          
     598 
    324599        for j from i >= j >= 0: 
    325600            counts[degrees[j]] -= 1 
    326601            output[counts[degrees[j]]] = self.wd_ents[start+j] 
    327          
     602 
    328603        max_location = counts[max_location] + start 
    329          
     604 
    330605        for j from 0 <= j <= i: 
    331606            self.wd_ents[start+j] = output[j] 
    332          
     607 
    333608        j = 1 
    334609        while j < self.nwords and counts[j] <= i: 
    … …  
    337612            self.wd_percolate(start + counts[j-1], start + counts[j] - 1) 
    338613            j += 1 
    339          
     614 
    340615        return max_location 
    341          
    342     cdef int refine(self, int k, int *col_alpha, int *wd_alpha, BinaryCodeGraph CG): 
     616 
     617    cdef int refine(self, int k, int *col_alpha, int *wd_alpha, BinaryCode CG): 
    343618        cdef int m = 0, j 
    344619        cdef int i, q, r, s, t 
    … …  
    416691            m += 1 
    417692        return invariant 
    418          
     693 
    419694    cdef void get_permutation(self, PartitionStack zeta, int *wd_gamma, int *col_gamma): 
    420695        cdef int i 
    … …  
    423698        for i from 0 <= i < self.nwords: 
    424699            wd_gamma[zeta.wd_ents[i]] = self.wd_ents[i] 
    425      
    426     cdef int cmp(self, PartitionStack other, BinaryCodeGraph CG): 
     700 
     701    cdef int cmp(self, PartitionStack other, BinaryCode CG): 
    427702        # if CG(self) > G(other): return 1 
    428703        # if CG(self) < G(other): return -1 
    … …  
    431706        for i from 0 <= i < CG.nwords: 
    432707            for j from 0 <= j < CG.ncols: 
    433                 k = CG.has_edge_bip(self.wd_ents[i], self.col_ents[j]) 
    434                 l = CG.has_edge_bip(other.wd_ents[i], other.col_ents[j]) 
     708                k = CG.is_one(self.wd_ents[i], self.col_ents[j]) 
     709                l = CG.is_one(other.wd_ents[i], other.col_ents[j]) 
    435710                if k ^ l: 
    436711                    return k - l 
    437712        return 0 
    438          
    439          
     713 
     714def classify(BinaryCode C, lab=True, verbosity=0): 
     715    """ 
     716    """ 
     717    cdef int i, j # local variables 
     718    cdef OrbitPartition Theta # keeps track of which vertices have been 
     719                              # discovered to be equivalent 
     720    cdef int index = 0, size = 1 
     721    cdef int L = 100 
     722    cdef int **Phi, **Omega 
     723    cdef int l = -1 
     724    cdef PartitionStack nu, zeta, rho 
     725    cdef int k_rho 
     726    cdef int k = 0 
     727    cdef int h = -1 
     728    cdef int hb 
     729    cdef int hh = 1 
     730    cdef int ht 
     731    cdef mpz_t *Lambda_mpz, *zf_mpz, *zb_mpz 
     732    cdef int hzf 
     733    cdef int hzb = -1 
     734    cdef unsigned int *basis_gamma 
     735    cdef int *col_gamma 
     736    cdef int *alpha 
     737    cdef int *v 
     738    cdef int *e 
     739    cdef int state 
     740    cdef int tvc, tvh, nwords = C.nwords, ncols = C.ncols, n = nwords + ncols, nrows = C.nrows 
     741 
     742    # trivial case 
     743    if ncols == 0: 
     744        return [], {} 
     745    elif nwords == 0 and ncols == 1: 
     746        return [], {0:0} 
     747    elif nwords == 0: 
     748        output1 = [] 
     749        dd = {} 
     750        for i from 0 <= i < ncols-1: 
     751            dd[i] = i 
     752            perm = range(ncols) 
     753            perm[i] = i+1 
     754            perm[i+1] = i 
     755            output1.append(perm) 
     756        dd[ncols-1] = ncols-1 
     757        return output1, dd 
     758 
     759    # allocate int pointers 
     760    Phi = <int **> sage_malloc(L * sizeof(int *)) 
     761    Omega = <int **> sage_malloc(L * sizeof(int *)) 
     762 
     763    # allocate GMP int pointers 
     764    Lambda_mpz = <mpz_t *> sage_malloc((ncols+2)*sizeof(mpz_t)) 
     765    zf_mpz     = <mpz_t *> sage_malloc((ncols+2)*sizeof(mpz_t)) 
     766    zb_mpz     = <mpz_t *> sage_malloc((ncols+2)*sizeof(mpz_t)) 
     767 
     768    # check for memory errors 
     769    if not (Phi and Omega and Lambda_mpz and zf_mpz and zb_mpz): 
     770        if Lambda_mpz: sage_free(Lambda_mpz) 
     771        if zf_mpz: sage_free(zf_mpz) 
     772        if zb_mpz: sage_free(zb_mpz) 
     773        if Phi: sage_free(Phi) 
     774        if Omega: sage_free(Omega) 
     775        raise MemoryError("Error allocating memory.") 
     776 
     777    # allocate int arrays 
     778    basis_gamma = <unsigned int *> sage_malloc(nrows*sizeof(unsigned int)) 
     779    col_gamma = <int *> sage_malloc(ncols*sizeof(int)) 
     780    Phi[0] = <int *> sage_malloc(L*ncols*sizeof(int)) 
     781    Omega[0] = <int *> sage_malloc(L*ncols*sizeof(int)) 
     782    alpha = <int *> sage_malloc(4*ncols*sizeof(int)) 
     783    v = <int *> sage_malloc(ncols*sizeof(int)) 
     784    e = <int *> sage_malloc(ncols*sizeof(int)) 
     785 
     786    # check for memory errors 
     787    if not (basis_gamma and col_gamma and Phi[0] and Omega[0] and alpha and v and e): 
     788        if basis_gamma: sage_free(basis_gamma) 
     789        if col_gamma: sage_free(col_gamma) 
     790        if Phi[0]: sage_free(Phi[0]) 
     791        if Omega[0]: sage_free(Omega[0]) 
     792        if alpha: sage_free(alpha) 
     793        if v: sage_free(v) 
     794        if e: sage_free(e) 
     795        sage_free(Lambda_mpz) 
     796        sage_free(zf_mpz) 
     797        sage_free(zb_mpz) 
     798        sage_free(Phi) 
     799        sage_free(Omega) 
     800        raise MemoryError("Error allocating memory.") 
     801 
     802    # setup double index arrays 
Note: See TracChangeset for help on using the changeset viewer.