Changeset 7762:e08b43bb53cf


Ignore:
Timestamp:
11/14/07 00:18:39 (6 years ago)
Author:
Robert L Miller <rlm@…>
Branch:
default
Message:

main loop builds

Location:
sage
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • sage/coding/binary_code.pxd

    r7761 r7762  
    1616 
    1717cdef class OrbitPartition: 
    18     cdef int nwords 
    1918    cdef int ncols 
    20     cdef int *wd_parent 
    21     cdef int *wd_rank 
    22     cdef int *wd_min_cell_rep 
    23     cdef int *wd_size 
    2419    cdef int *col_parent 
    2520    cdef int *col_rank 
    … …  
    2722    cdef int *col_size 
    2823 
    29     cdef int wd_find(self, int) 
    30     cdef void wd_union(self, int, int) 
    3124    cdef int col_find(self, int) 
    3225    cdef void col_union(self, int, int) 
    33     cdef int merge_perm(self, int *, int *) 
     26    cdef int merge_perm(self, int *) 
    3427 
    3528cdef class PartitionStack: 
    … …  
    4437    cdef int radix 
    4538    cdef int flag 
     39    cdef int v 
    4640    cdef int *col_degs   # 
    4741    cdef int *col_counts # 
    … …  
    7266cdef class BinaryCodeClassifier: 
    7367    cdef int *ham_wts 
     68    cdef int L 
     69    cdef int *Phi 
     70    cdef int *Omega 
     71    cdef int *W 
     72    cdef int radix 
     73    cdef int *Lambda1, *Lambda2, *Lambda3 
     74    cdef int *b_gamma, *c_gamma 
     75    cdef int *alpha 
     76    cdef int alpha_size 
     77    cdef int *v, *e 
     78    cdef int *aut_gp_gens, *labeling 
     79    cdef int aut_gp_index, aut_gens_size 
     80 
     81    cdef void record_automorphism(self, int *, int) 
     82    cdef void aut_gp_and_can_label(self, BinaryCode, int) 
    7483 
    7584 
    … …  
    8089 
    8190 
    82  

  • sage/coding/binary_code.pyx

    r7761 r7762  
    399399 
    400400    """ 
    401     def __new__(self, nrows, ncols): 
    402         cdef int nwords, word 
     401    def __new__(self, ncols): 
    403402        cdef int col 
    404         nwords = (1 << nrows) 
    405         self.nwords = nwords 
    406403        self.ncols = ncols 
    407         self.wd_parent =       <int *> sage_malloc( nwords * sizeof(int) ) 
    408         self.wd_rank =         <int *> sage_malloc( nwords * sizeof(int) ) 
    409         self.wd_min_cell_rep = <int *> sage_malloc( nwords * sizeof(int) ) 
    410         self.wd_size =         <int *> sage_malloc( nwords * sizeof(int) ) 
    411404        self.col_parent =       <int *> sage_malloc( ncols * sizeof(int) ) 
    412405        self.col_rank =         <int *> sage_malloc( ncols * sizeof(int) ) 
    413406        self.col_min_cell_rep = <int *> sage_malloc( ncols * sizeof(int) ) 
    414407        self.col_size =         <int *> sage_malloc( ncols * sizeof(int) ) 
    415         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): 
    416             if self.wd_parent: sage_free(self.wd_parent) 
    417             if self.wd_rank: sage_free(self.wd_rank) 
    418             if self.wd_min_cell_rep: sage_free(self.wd_min_cell_rep) 
    419             if self.wd_size: sage_free(self.wd_size) 
    420             if self.col_parent: sage_free(self.col_parent) 
    421             if self.col_rank: sage_free(self.col_rank) 
     408        if not (self.col_parent and self.col_rank and self.col_min_cell_rep and self.col_size): 
     409            if self.col_parent:       sage_free(self.col_parent) 
     410            if self.col_rank:         sage_free(self.col_rank) 
    422411            if self.col_min_cell_rep: sage_free(self.col_min_cell_rep) 
    423             if self.col_size: sage_free(self.col_size) 
     412            if self.col_size:         sage_free(self.col_size) 
    424413            raise MemoryError("Memory.") 
    425         for word from 0 <= word < nwords: 
    426             self.wd_parent[word] = word 
    427             self.wd_rank[word] = 0 
    428             self.wd_min_cell_rep[word] = word 
    429             self.wd_size[word] = 1 
    430414        for col from 0 <= col < ncols: 
    431415            self.col_parent[col] = col 
    … …  
    435419 
    436420    def __dealloc__(self): 
    437         sage_free(self.wd_parent) 
    438         sage_free(self.wd_rank) 
    439         sage_free(self.wd_min_cell_rep) 
    440         sage_free(self.wd_size) 
    441421        sage_free(self.col_parent) 
    442422        sage_free(self.col_rank) 
    … …  
    462442        cdef int i 
    463443        cdef int j 
    464         s = 'OrbitPartition on %d words and %d columns. Data:\nWords:\n'%(self.nwords, self.ncols) 
    465         for i from 0 <= i < self.nwords: 
    466             s += '%d,'%self.wd_parent[i] 
     444        s = 'OrbitPartition on %d columns. Data:\n'%(self.ncols) 
    467445        s = s[:-1] + '\nColumns:\n' 
    468446        for j from 0 <= j < self.ncols: 
    469447            s += '%d,'%self.col_parent[j] 
    470448        return s[:-1] 
    471  
    472     def _wd_find(self, word): 
    473         """ 
    474         Returns the root of word. 
    475  
    476         EXAMPLE: 
    477             sage: import sage.coding.binary_code 
    478             sage: from sage.coding.binary_code import * 
    479             sage: O = OrbitPartition(4, 8) 
    480             sage: O 
    481             OrbitPartition on 16 words and 8 columns. Data: 
    482             Words: 
    483             0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
    484             Columns: 
    485             0,1,2,3,4,5,6,7 
    486             sage: O._wd_find(12) 
    487             12 
    488  
    489         """ 
    490         return self.wd_find(word) 
    491  
    492     cdef int wd_find(self, int word): 
    493         if self.wd_parent[word] == word: 
    494             return word 
    495         else: 
    496             self.wd_parent[word] = self.wd_find(self.wd_parent[word]) 
    497             return self.wd_parent[word] 
    498  
    499     def _wd_union(self, x, y): 
    500         """ 
    501         Join the cells containing x and y. 
    502  
    503         EXAMPLE: 
    504             sage: import sage.coding.binary_code 
    505             sage: from sage.coding.binary_code import * 
    506             sage: O = OrbitPartition(4, 8) 
    507             sage: O 
    508             OrbitPartition on 16 words and 8 columns. Data: 
    509             Words: 
    510             0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 
    511             Columns: 
    512             0,1,2,3,4,5,6,7 
    513             sage: O._wd_union(1,10) 
    514             sage: O 
    515             OrbitPartition on 16 words and 8 columns. Data: 
    516             Words: 
    517             0,1,2,3,4,5,6,7,8,9,1,11,12,13,14,15 
    518             Columns: 
    519             0,1,2,3,4,5,6,7 
    520             sage: O._wd_find(10) 
    521             1 
    522  
    523         """ 
    524         self.wd_union(x, y) 
    525  
    526     cdef void wd_union(self, int x, int y): 
    527         cdef int x_root, y_root 
    528         x_root = self.wd_find(x) 
    529         y_root = self.wd_find(y) 
    530         if self.wd_rank[x_root] > self.wd_rank[y_root]: 
    531             self.wd_parent[y_root] = x_root 
    532             self.wd_min_cell_rep[y_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root]) 
    533             self.wd_size[y_root] += self.wd_size[x_root] 
    534         elif self.wd_rank[x_root] < self.wd_rank[y_root]: 
    535             self.wd_parent[x_root] = y_root 
    536             self.wd_min_cell_rep[x_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root]) 
    537             self.wd_size[x_root] += self.wd_size[y_root] 
    538         elif x_root != y_root: 
    539             self.wd_parent[y_root] = x_root 
    540             self.wd_min_cell_rep[y_root] = min(self.wd_min_cell_rep[x_root],self.wd_min_cell_rep[y_root]) 
    541             self.wd_size[y_root] += self.wd_size[x_root] 
    542             self.wd_rank[x_root] += 1 
    543449 
    544450    def _col_find(self, col): 
    … …  
    614520            self.col_rank[x_root] += 1 
    615521 
    616     def _merge_perm(self, col_gamma, wd_gamma): 
     522    def _merge_perm(self, col_gamma): 
    617523        """ 
    618524        Merges the cells of self under the given permutation. If gamma[a] = b, 
    … …  
    642548        cdef int i 
    643549        cdef int *_col_gamma 
    644         cdef int *_wd_gamma 
    645         _wd_gamma = <int *> sage_malloc(self.nwords * sizeof(int)) 
    646550        _col_gamma = <int *> sage_malloc(self.ncols * sizeof(int)) 
    647         if not (_col_gamma and _wd_gamma): 
    648             if _wd_gamma: sage_free(_wd_gamma) 
    649             if _col_gamma: sage_free(_col_gamma) 
     551        if not _col_gamma: 
    650552            raise MemoryError("Memory.") 
    651         for i from 0 <= i < self.nwords: 
    652             _wd_gamma[i] = wd_gamma[i] 
    653553        for i from 0 <= i < self.ncols: 
    654554            _col_gamma[i] = col_gamma[i] 
    655         result = self.merge_perm(_col_gamma, _wd_gamma) 
     555        result = self.merge_perm(_col_gamma) 
    656556        sage_free(_col_gamma) 
    657         sage_free(_wd_gamma) 
    658557        return result 
    659558 
    660     cdef int merge_perm(self, int *col_gamma, int *wd_gamma): 
     559    cdef int merge_perm(self, int *col_gamma): 
    661560        cdef int i, gamma_i_root 
    662561        cdef int j, gamma_j_root, return_value = 0 
    663         cdef int *self_wd_parent = self.wd_parent 
    664562        cdef int *self_col_parent = self.col_parent 
    665         for i from 0 <= i < self.nwords: 
    666             if self_wd_parent[i] == i: 
    667                 gamma_i_root = self.wd_find(wd_gamma[i]) 
    668                 if gamma_i_root != i: 
    669                     return_value = 1 
    670                     self.wd_union(i, gamma_i_root) 
    671563        for j from 0 <= j < self.ncols: 
    672564            if self_col_parent[j] == j: 
    … …  
    697589            self.nwords = other.nwords 
    698590            self.ncols = other.ncols 
     591 
    699592        self.radix = 8*sizeof(int) 
    700593        self.flag = (1 << (self.radix-1)) 
    … …  
    717610            and self.col_degs and self.col_counts and self.col_output \ 
    718611            and self.wd_degs  and self.wd_counts  and self.wd_output): 
    719             if self.wd_ents:    sage_free(self.wd_ents) 
    720             if self.wd_lvls:    sage_free(self.wd_lvls) 
    721             if self.col_ents:   sage_free(self.col_ents) 
    722             if self.col_lvls:   sage_free(self.col_lvls) 
    723             if self.col_degs:   sage_free(self.col_degs) 
    724             if self.col_counts: sage_free(self.col_counts) 
    725             if self.col_output: sage_free(self.col_output) 
    726             if self.wd_degs:    sage_free(self.wd_degs) 
    727             if self.wd_counts:  sage_free(self.wd_counts) 
    728             if self.wd_output:  sage_free(self.wd_output) 
     612            if self.wd_ents:         sage_free(self.wd_ents) 
     613            if self.wd_lvls:         sage_free(self.wd_lvls) 
     614            if self.col_ents:        sage_free(self.col_ents) 
     615            if self.col_lvls:        sage_free(self.col_lvls) 
     616            if self.col_degs:        sage_free(self.col_degs) 
     617            if self.col_counts:      sage_free(self.col_counts) 
     618            if self.col_output:      sage_free(self.col_output) 
     619            if self.wd_degs:         sage_free(self.wd_degs) 
     620            if self.wd_counts:       sage_free(self.wd_counts) 
     621            if self.wd_output:       sage_free(self.wd_output) 
    729622            raise MemoryError("Memory.") 
    730623 
    … …  
    11821075        # nontrivial cell 
    11831076        j = location 
     1077        self.v = j 
    11841078        while True: 
    11851079            if self_col_lvls[j] <= k: break 
    … …  
    13281222        """ 
    13291223        Split column v out, placing it before the rest of the cell it was in. 
     1224        Returns the location of the split column. 
    13301225 
    13311226        EXAMPLE: 
    … …  
    17891684                    else: j = i 
    17901685            m += 1 
    1791         return invariant 
     1686        if invariant != -1: 
     1687            return invariant 
     1688        else: 
     1689            return 0 
    17921690 
    17931691    def _clear(self, k): 
    … …  
    20941992            col_gamma[other_col_ents[i]] = self_col_ents[i] 
    20951993 
    2096 ################################################################################ 
    2097 ################################################################################ 
    2098 ################################################################################ 
    2099  
    21001994cdef class BinaryCodeClassifier: 
    21011995 
    21021996    def __new__(self): 
     1997        self.radix = sizeof(int) << 3 
    21031998        self.ham_wts = hamming_weights() 
     1999        self.L = 100 # memory limit for Phi and Omega 
     2000        self.alpha_size = 65536 + self.radix 
     2001        self.Phi =         <int *> sage_malloc( self.L           * sizeof(int) ) 
     2002        self.Omega =       <int *> sage_malloc( self.L           * sizeof(int) ) 
     2003        self.W =           <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2004        self.Lambda1 =     <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2005        self.Lambda2 =     <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2006        self.Lambda3 =     <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2007        self.b_gamma =     <int *> sage_malloc( (self.radix-1)   * sizeof(int) ) 
     2008        self.c_gamma =     <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2009        self.alpha =       <int *> sage_malloc( self.alpha_size  * sizeof(int) ) 
     2010        self.v =           <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2011        self.e =           <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2012        self.aut_gp_gens = <int *> sage_malloc( self.radix * 100 * sizeof(int) ) 
     2013        self.aut_gens_size = self.radix * 100 
     2014        self.labeling =    <int *> sage_malloc( self.radix       * sizeof(int) ) 
     2015        if not (self.Phi and self.Omega and self.W and self.Lambda1 and self.Lambda2 and self.Lambda3 \ 
     2016            and self.b_gamma and self.c_gamma and self.alpha and self.v and self.e and self.aut_gp_gens \ 
     2017            and self.labeling): 
     2018            if self.Phi:          sage_free(self.Phi) 
     2019            if self.Omega:        sage_free(self.Omega) 
     2020            if self.W:            sage_free(self.W) 
     2021            if self.Lambda1:      sage_free(self.Lambda1) 
     2022            if self.Lambda2:      sage_free(self.Lambda2) 
     2023            if self.Lambda3:      sage_free(self.Lambda3) 
     2024            if self.b_gamma:      sage_free(self.b_gamma) 
     2025            if self.c_gamma:      sage_free(self.c_gamma) 
     2026            if self.alpha:        sage_free(self.alpha) 
     2027            if self.v:            sage_free(self.v) 
     2028            if self.e:            sage_free(self.e) 
     2029            if self.aut_gp_gens:  sage_free(self.aut_gp_gens) 
     2030            if self.labeling:     sage_free(self.labeling) 
     2031            raise MemoryError("Memory.") 
    21042032 
    21052033    def __dealloc__(self): 
    21062034        sage_free(self.ham_wts) 
    2107  
    2108  
    2109  
    2110  
    2111  
    2112  
    2113  
    2114  
    2115  
    2116  
    2117  
    2118  
    2119  
    2120  
    2121  
    2122  
    2123  
    2124  
    2125  
    2126 #def classify(BinaryCode C, lab=True, verbosity=0): 
    2127 #    """ 
    2128 #    """ 
    2129 #    cdef int i, j # local variables 
    2130 #    cdef OrbitPartition Theta # keeps track of which vertices have been 
    2131 #                              # discovered to be equivalent 
    2132 #    cdef int index = 0, size = 1 
    2133 #    cdef int L = 100 
    2134 #    cdef int **Phi, **Omega 
    2135 #    cdef int l = -1 
    2136 #    cdef PartitionStack nu, zeta, rho 
    2137 #    cdef int k_rho 
    2138 #    cdef int k = 0 
    2139 #    cdef int h = -1 
    2140 #    cdef int hb 
    2141 #    cdef int hh = 1 
    2142 #    cdef int ht 
    2143 #    cdef mpz_t *Lambda_mpz, *zf_mpz, *zb_mpz 
    2144 #    cdef int hzf 
    2145 #    cdef int hzb = -1 
    2146 #    cdef int *basis_gamma 
    2147 #    cdef int *col_gamma 
    2148 #    cdef int *alpha 
    2149 #    cdef int *v 
    2150 #    cdef int *e 
    2151 #    cdef int state 
    2152 #    cdef int tvc, tvh, nwords = C.nwords, ncols = C.ncols, n = nwords + ncols, nrows = C.nrows 
    2153  
    2154 #    # trivial case 
    2155 #    if ncols == 0: 
    2156 #        return [], {} 
    2157 #    elif nwords == 0 and ncols == 1: 
    2158 #        return [], {0:0} 
    2159 #    elif nwords == 0: 
    2160 #        output1 = [] 
    2161 #        dd = {} 
    2162 #        for i from 0 <= i < ncols-1: 
    2163 #            dd[i] = i 
    2164 #            perm = range(ncols) 
    2165 #            perm[i] = i+1 
    2166 #            perm[i+1] = i 
    2167 #            output1.append(perm) 
    2168 #        dd[ncols-1] = ncols-1 
    2169 #        return output1, dd 
    2170  
    2171 #    # allocate int pointers 
    2172 #    Phi = <int **> sage_malloc(L * sizeof(int *)) 
    2173 #    Omega = <int **> sage_malloc(L * sizeof(int *)) 
    2174  
    2175 #    # allocate GMP int pointers 
    2176 #    Lambda_mpz = <mpz_t *> sage_malloc((ncols+2)*sizeof(mpz_t)) 
    2177 #    zf_mpz     = <mpz_t *> sage_malloc((ncols+2)*sizeof(mpz_t)) 
    2178 #    zb_mpz     = <mpz_t *> sage_malloc((ncols+2)*sizeof(mpz_t)) 
    2179  
    2180 #    # check for memory errors 
    2181 #    if not (Phi and Omega and Lambda_mpz and zf_mpz and zb_mpz): 
    2182 #        if Lambda_mpz: sage_free(Lambda_mpz) 
    2183 #        if zf_mpz: sage_free(zf_mpz) 
    2184 #        if zb_mpz: sage_free(zb_mpz) 
    2185 #        if Phi: sage_free(Phi) 
    2186 #        if Omega: sage_free(Omega) 
    2187 #        raise MemoryError("Error allocating memory.") 
    2188  
    2189 #    # allocate int arrays 
    2190 #    basis_gamma = <int *> sage_malloc(nrows*sizeof(int)) 
    2191 #    col_gamma = <int *> sage_malloc(ncols*sizeof(int)) 
    2192 #    Phi[0] = <int *> sage_malloc(L*ncols*sizeof(int)) 
    2193 #    Omega[0] = <int *> sage_malloc(L*ncols*sizeof(int)) 
    2194 #    alpha = <int *> sage_malloc(4*ncols*sizeof(int)) 
    2195 #    v = <int *> sage_malloc(ncols*sizeof(int)) 
    2196 #    e = <int *> sage_malloc(ncols*sizeof(int)) 
    2197  
    2198 #    # check for memory errors 
    2199 #    if not (basis_gamma and col_gamma and Phi[0] and Omega[0] and alpha and v and e): 
    2200 #        if basis_gamma: sage_free(basis_gamma) 
    2201 #        if col_gamma: sage_free(col_gamma) 
    2202 #        if Phi[0]: sage_free(Phi[0]) 
    2203 #        if Omega[0]: sage_free(Omega[0]) 
    2204 #        if alpha: sage_free(alpha) 
    2205 #        if v: sage_free(v) 
    2206 #        if e: sage_free(e) 
    2207 #        sage_free(Lambda_mpz) 
    2208 #        sage_free(zf_mpz) 
    2209 #        sage_free(zb_mpz) 
    2210 #        sage_free(Phi) 
    2211 #        sage_free(Omega) 
    2212 #        raise MemoryError("Error allocating memory.") 
    2213  
    2214 #    # setup double index arrays 
    2215  
    2216  
    2217  
    2218  
    2219  
    2220  
    2221  
    2222  
    2223  
     2035        sage_free(self.Phi) 
     2036        sage_free(self.Omega) 
     2037        sage_free(self.W) 
     2038        sage_free(self.Lambda1) 
     2039        sage_free(self.Lambda2) 
     2040        sage_free(self.Lambda3) 
     2041        sage_free(self.b_gamma) 
     2042        sage_free(self.c_gamma) 
     2043        sage_free(self.alpha) 
     2044        sage_free(self.v) 
     2045        sage_free(self.e) 
     2046        sage_free(self.aut_gp_gens) 
     2047        sage_free(self.labeling) 
     2048 
     2049    cdef void record_automorphism(self, int *gamma, int ncols): 
     2050        cdef int i, j 
     2051        if self.aut_gp_index + ncols > self.aut_gens_size: 
     2052            self.aut_gens_size *= 2 
     2053            self.aut_gp_gens = <int *> sage_realloc( self.aut_gp_gens, self.aut_gens_size ) 
     2054            if not self.aut_gp_gens: 
     2055                raise MemoryError("Memory.") 
     2056        j = self.aut_gp_index 
     2057        for i from 0 <= i < ncols: 
     2058            self.aut_gp_gens[i+j] = gamma[i] 
     2059        self.aut_gp_index += ncols 
     2060 
     2061    def _aut_gp_and_can_label(self, C, verbosity=0): 
     2062        cdef int i, j 
     2063        self.aut_gp_and_can_label(C, verbosity) 
     2064        i = 0 
     2065        py_aut_gp_gens = [] 
     2066        while self.aut_gp_gens[i] != -1: 
     2067            gen = [self.aut_gp_gens[i+j] for j from 0 <= j < C.ncols] 
     2068            py_aut_gp_gens.append(gen) 
     2069            i += C.ncols 
     2070        py_labeling = [self.labeling[i] for i from 0 <= i < C.ncols] 
     2071        return py_aut_gp_gens, py_labeling 
     2072 
     2073    cdef void aut_gp_and_can_label(self, BinaryCode C, int verbosity): 
     2074 
     2075        # allocate/set up variables: 
     2076        cdef int i, j, ii, jj # local variables 
     2077        cdef OrbitPartition Theta # keeps track of which vertices have been 
     2078                                  # discovered to be equivalent 
     2079        cdef int index = 0, size = 1    # Define $\Gamma^{(-1)} := \text{Aut}(C)$, and 
     2080                                        # $\Gamma^{(i)} := \Gamma^{(-1)}_{v_0,...,v_i}$. 
     2081                                        # Then index = $|\Gamma^{(k-1)}|/|\Gamma^{(k)}|$ at (POINT A) 
     2082                                        # and size = $|\Gamma^{(k-1)}|$ at (POINT A) and (POINT B). 
     2083        cdef int *Phi   = self.Phi      # Phi stores the fixed point sets of each automorphism 
     2084        cdef int *Omega = self.Omega    # Omega stores the minimal elements of each cell of the orbit partition 
     2085        cdef int l = -1 # current index for storing values in Phi and Omega- we start at -1 so that when 
     2086                        # we increment first, the first place we write to is 0. 
     2087        cdef int *W = self.W    # for each k, W[k] is a list (as int mask) of the vertices to be searched down from 
     2088                                # the current partition, at k. Phi and Omega are ultimately used to make the size of 
     2089                                # W as small as possible 
     2090        cdef PartitionStack nu, zeta, rho 
     2091        cdef int k_rho  # the number of partitions in rho 
     2092        cdef int k = 0  # the number of partitions in nu 
     2093        cdef int h = -1 # longest common ancestor of zeta and nu: 
     2094                        # zeta[h] == nu[h], zeta[h+1] != nu[h+1] 
     2095                        # -1 indicates that zeta is not yet defined 
     2096        cdef int hb     # longest common ancestor of rho and nu: 
     2097                        # rho[hb] == nu[hb], rho[hb+1] != nu[hb+1] 
     2098        cdef int hh = 1 # the height of the oldest ancestor of nu 
     2099                        # satisfying Lemma 2.25 in [1]: 
     2100                        # if nu does not satisfy it at k, then hh = k 
     2101        cdef int ht # smallest such that all descendants of zeta[ht] are equivalent 
     2102        cdef int *Lambda = self.Lambda1             # for tracking indicator values- zf and zb are 
     2103        cdef int *zf__Lambda_zeta = self.Lambda2    # indicator vectors remembering Lambda[k] for 
     2104        cdef int *zb__Lambda_rho = self.Lambda3     # zeta and rho, respectively 
     2105        cdef int qzb    # keeps track of Lambda[k] ><= zb[k] 
     2106        cdef int hzf__h_zeta      # the max height for which Lambda and zf agree 
     2107        cdef int hzb__h_rho = -1  # the max height for which Lambda and zb agree 
     2108        cdef int *basis_gamma = self.b_gamma, *col_gamma = self.c_gamma # used for storing permutations 
     2109        cdef int tvc, tvh, nwords = C.nwords, ncols = C.ncols, nrows = C.nrows 
     2110        # TODO: Explain the variables tvc, tvh 
     2111        cdef int *alpha # for storing pointers to cells of nu[k] 
     2112        cdef int *v = self.v    # list of vertices determining nu 
     2113        cdef int *e = self.e    # 0 or 1, see states 12 and 17 
     2114        cdef int state = 1  # keeps track of position in algorithm - see sage/graphs/graph_isom.pyx, search for "STATE DIAGRAM" 
     2115        cdef int *ham_wts = self.ham_wts 
     2116        self.aut_gp_index = 0 
     2117        while self.alpha_size < nrows + ncols: 
     2118            self.alpha_size *= 2 
     2119            self.alpha = <int *> sage_realloc(self.alpha, self.alpha_size * sizeof(int) ) 
     2120            if not self.alpha: 
     2121                raise MemoryError("Memory.") 
     2122        alpha = self.alpha # think of alpha as of length exactly nrows + ncols 
     2123        nu =    PartitionStack(nrows, ncols) 
     2124        Theta = OrbitPartition(ncols) 
     2125 
     2126        # trivial case 
     2127        if ncols == 0 or nrows == 0: 
     2128            raise NotImplementedError("Must supply a nontrivial code.") 
     2129 
     2130        while state != -1: 
     2131 
     2132            if state == 1: # Entry point: once only 
     2133                alpha[0] = 0 
     2134                alpha[1] = nu.flag 
     2135                nu.refine(k, alpha, 2, C, ham_wts) 
     2136                if nu.sat_225(k): hh = k 
     2137                if nu.is_discrete(k): state = 18; continue 
     2138 
     2139                # store the first smallest nontrivial cell in W[k], and set v[k] 
     2140                # equal to its minimum element 
     2141                W[k] = nu.first_smallest_nontrivial(k) 
     2142                v[k] = nu.v # stored during first_smallest_nontrivial 
     2143                Lambda[k] = 0 
     2144                e[k] = 0 
     2145                state = 2 
     2146 
     2147            elif state == 2: # Move down the search tree one level by refining nu: 
     2148                             # split out a column, and refine nu against it 
     2149                k += 1 
     2150                nu.clear(k) 
     2151                alpha[0] = nu.split_column(v[k-1], k) 
     2152                Lambda[k] = nu.refine(k, alpha, 1, C, ham_wts) # store the invariant to Lambda[k] 
     2153                # only if this is the first time moving down the search tree: 
     2154                if h == -1: state = 5; continue 
     2155 
     2156                # update hzf__h_zeta 
     2157                if hzf__h_zeta == k-1 and Lambda[k] == zf__Lambda_zeta[k]: hzf__h_zeta = k 
     2158                # update qzb 
     2159                if zb__Lambda_rho[k] == -1 or Lambda[k] < zb__Lambda_rho[k]: 
     2160                    qzb = -1 
     2161                elif Lambda[k] > zb__Lambda_rho[k]: 
     2162                    qzb = 1 
     2163                else: 
     2164                    qzb = 0 
     2165                # update hzb 
     2166                if hzb__h_rho == k-1 and qzb == 0: hzb__h_rho = k 
     2167                # if Lambda[k] > zb[k], then zb[k] := Lambda[k] 
     2168                # (zb keeps track of the indicator invariants corresponding to 
     2169                # rho, the closest canonical leaf so far seen- if Lambda is 
     2170                # bigger, then rho must be about to change 
     2171                if qzb > 0: zb__Lambda_rho[k] = Lambda[k] 
     2172                state = 3 
     2173 
     2174            elif state == 3: # attempt to rule out automorphisms while moving down the tree 
     2175                # if k > hzf, then we know that nu currently does not look like zeta, the first 
     2176                # terminal node encountered, thus there is no automorphism to discover. If qzb < 0, 
     2177                # i.e. Lambda[k] < zb[k], then the indicator is not maximal, and we can't reach a 
     2178                # canonical leaf. If neither of these is the case, then proceed to state 6. 
     2179                if hzf__h_zeta <= k or qzb >= 0: state = 4 
     2180                else: state = 6 
     2181 
     2182            elif state == 4: # at this point we have -not- ruled out the presence of automorphisms 
     2183                if nu.is_discrete(k): state = 7; continue # we have a terminal node, so process it 
     2184 
     2185                # otherwise, prepare to split out another column: 
     2186                # store the first smallest nontrivial cell in W[k], and set v[k] 
     2187                # equal to its minimum element 
     2188                W[k] = nu.first_smallest_nontrivial(k) 
     2189                v[k] = nu.v # stored during first_smallest_nontrivial 
     2190                if not nu.sat_225(k): hh = k + 1 
     2191                e[k] = 0 # see state 12 and 17 
     2192                state = 2 # continue down the tree 
     2193 
     2194            elif state == 5: # same as state 3, but in the case where we haven't yet defined zeta 
     2195                             # i.e. this is our first time down the tree. Once we get to the bottom, 
     2196                             # we will have zeta = nu = rho, so we do: 
     2197                zf__Lambda_zeta[k] = Lambda[k] 
     2198                zb__Lambda_rho[k] = Lambda[k] 
     2199                state = 4 
     2200 
     2201            elif state == 6: # at this stage, there is no reason to continue downward, so backtrack 
     2202                j = k 
     2203                # return to the longest ancestor nu[i] of nu that could have a 
     2204                # descendant equivalent to zeta or could improve on rho. 
     2205                # All terminal nodes descending from nu[hh] are known to be 
     2206                # equivalent, so i < hh. Also, if i > hzb, none of the 
     2207                # descendants of nu[i] can improve rho, since the indicator is 
     2208                # off (Lambda(nu) < Lambda(rho)). If i >= ht, then no descendant 
     2209                # of nu[i] is equivalent to zeta (see [1, p67]). 
     2210                if ht-1 > hzb__h_rho: 
     2211                    if ht-1 < hh-1: 
     2212                        k = ht-1 
     2213                    else: 
     2214                        k = hh-1 
     2215                else: 
     2216                    if hzb__h_rho < hh-1: 
     2217                        k = hzb__h_rho 
     2218                    else: 
     2219                        k = hh-1 
     2220                # TODO: is the following line necessary? 
     2221                if k == -1: k = 0 
     2222                # if j == hh, then all nodes lower than our current position are equivalent, so bail out 
     2223                if j == hh: state = 13; continue 
     2224 
     2225                # recall hh: the height of the oldest ancestor of zeta for which Lemma 2.25 is 
     2226                # satsified, which implies that all terminal nodes descended from there are equivalent. 
     2227                # If we are looking at such a node, then the partition at nu[hh] can be used for later 
     2228                # pruning, so we store its fixed set and a set of representatives of its cells. 
     2229                if l < self.L-1: l += 1 
     2230                Omega[l] = nu.min_cell_reps(hh) 
     2231                Phi[l] = nu.fixed_cols(Omega[l], hh) 
     2232                state = 12 
     2233 
     2234            elif state == 7: # we have just arrived at a terminal node of the search tree T(G, Pi) 
     2235                # if this is the first terminal node, go directly to 18, to 
     2236                # process zeta 
     2237                if h == -1: state = 18; continue 
     2238 
     2239                # hzf is the extremal height of ancestors of both nu and zeta, so if k < hzf, nu is not 
     2240                # equivalent to zeta, i.e. there is no automorphism to discover. 
     2241                if k < hzf: state = 8; continue 
     2242 
     2243                nu.get_permutation(zeta, basis_gamma, col_gamma, ham_wts) 
     2244                # if C^gamma == C, the permutation is an automorphism, goto 10 
     2245                if C.is_automorphism(col_gamma, basis_gamma): 
     2246                    state = 10 
     2247                else: 
     2248                    state = 8 
     2249 
     2250            elif state == 8: # we have just ruled out the presence of automorphism and have not yet 
     2251                             # considered whether nu improves on rho 
     2252                # if qzb < 0, then rho already has larger indicator tuple 
     2253                if qzb < 0: state = 6; continue 
     2254 
     2255                # if Lambda[k] > zb[k] or nu is shorter than rho, then we have an improvement for rho 
     2256                if (qzb > 0) or (k < k_rho): state = 9; continue 
     2257 
     2258                # now Lambda[k] == zb[k] and k == k_rho, so we appeal to an enumeration: 
     2259                j = nu.cmp(rho, C) 
     2260                # if C(nu) > C(rho), we have a new label, goto 9 
     2261                if j > 0: state = 9; continue 
     2262 
     2263                # if C(nu) < C(rho), no new label, goto 6 
     2264                if j < 0: state = 6; continue 
     2265 
     2266                # if C(nu) == C(rho), get the automorphism and goto 10 
     2267                nu.find_basis(ham_wts) 
     2268                rho.get_permutation(nu, basis_gamma, col_gamma, ham_wts) 
     2269                state = 10 
     2270 
     2271            elif state == 9: # nu is a better guess at the canonical label than rho 
     2272                rho = PartitionStack(nu) 
     2273                k_rho = k 
     2274                qzb = 0 
     2275                hb = k 
     2276                hzb__h_rho = k 
     2277                # set zb[k+1] = Infinity 
     2278                zb__Lambda_rho[k+1] = -1 
     2279                state = 6 
     2280 
     2281            elif state == 10: # we have an automorphism to process 
     2282                # increment l 
     2283                if l < self.L-1: l += 1 
     2284                # store information about the automorphism to Omega and Phi 
     2285                # Omega stores the minimum cell representatives 
     2286                Omega[l] = ~0 
     2287                i = 0 
     2288                while i < ncols: 
     2289                    j = col_gamma[i]         # i is a minimum 
     2290                    while j != i:            # cell rep, 
     2291                        Omega[l] ^= (1<<j)   # so cancel 
     2292                        j = col_gamma[j]     # cellmates 
     2293                    while i < ncols and not Omega[l]&(1<<i): # find minimal element 
     2294                        i += 1                               # of next cell 
     2295                # Phi stores the columns fixed by the automorphism 
     2296                Phi[l] = 0 
     2297                for i from 0 <= i < ncols: 
     2298                    if col_gamma[i] == i: 
     2299                        Phi[l] ^= (1 << i) 
     2300                # Now incorporate the automorphism into Theta 
     2301                j = Theta.merge_perm(col_gamma) 
     2302                # j stores whether anything happened or not- if not, then the automorphism we have 
     2303                # discovered is already in the subgroup spanned by the generators we have output 
     2304                if not j: state = 11; continue 
     2305 
     2306                # otherwise, we have a new generator, so record it: 
     2307                self.record_automorphism(col_gamma, ncols) 
     2308 
     2309                # The variable tvc was set to be the minimum element of W[k] the last time the 
     2310                # algorithm came up to meet zeta. At this point, we were considering the new 
     2311                # possibilities for descending away from zeta at this level. 
     2312                if Theta.col_min_cell_rep[Theta.col_find(tvc)] == tvc: 
     2313                    # if this is still a minimum cell representative of Theta, even in light 
     2314                    # of this new automorphism, then the current branch off of zeta hasn't been 
     2315                    # found equivalent to one already searched yet, so there may still be a 
     2316                    # better canonical label downward. 
     2317                    state = 11; continue 
     2318 
     2319                # Otherwise, proceed to where zeta meets nu: 
     2320                k = h 
     2321                state = 13 
     2322 
     2323            elif state == 11: # We have just found a new automorphism, and deduced that there may 
     2324                # be a better canonical label below the current branch off of zeta. So go to where 
     2325                # nu meets rho 
     2326                k = hb 
     2327                state = 12 
     2328 
     2329            elif state == 12: # Coming here from either state 6 or 11, the algorithm has discovered 
     2330                              # some new information. 11 came from 10, where a new line in Omega and 
     2331                              # Phi was just recorded, and 6 stored information about implicit auto- 
     2332                              # morphisms in Omega and Phi 
     2333                if e[k] == 1: 
     2334                    # this means that the algorithm has come upward to this position (in state 17) 
     2335                    # before, so we have already intersected W[k] with the bulk of Omega and Phi, but 
     2336                    # we should still catch up with the latest ones 
     2337                    W[k] &= Omega[l] 
     2338                state = 13 
     2339 
     2340            elif state == 13: # hub state 
     2341                if k == -1: state = -1; continue # exit point 
     2342 
     2343                if k > h: state = 17; continue # we are still on the same principal branch from zeta 
     2344 
     2345                if k == h: state = 14; continue # update the stabilizer index and check for new splits, 
     2346                                                # since we have returned to a partition of zeta 
     2347                # otherwise k < h, hence we have just backtracked up zeta, and are one level closer to done 
     2348                h = k 
     2349                tvc = 0 
     2350                while not (1 << tvc) & W[k]: 
     2351                    tvc += 1 
     2352                # now tvc points to the minimal cell representative of W[k] 
     2353                state = 14 
     2354 
     2355            elif state == 14: # see if there are any more splits to make from this level of zeta (see state 17) 
     2356                if Theta.col_find(v[k]) == Theta.col_find(tvc): 
     2357                    index += 1 
     2358                    # keep tabs on how many elements are in the same cell of Theta as tvc 
     2359                # find the next split 
     2360                i = v[k] + 1 
     2361                while i < ncols and not (1 << i) & W[k]: 
     2362                    i += 1 
     2363                if i < ncols: 
     2364                    v[k] = i 
     2365                else: 
     2366                    # there is no new split at this level 
     2367                    v[k] = -1 
     2368                    state = 16; continue 
     2369 
     2370                # new split column better be a minimal representative in Theta, or wasted effort 
     2371                if Theta.col_min_cell_rep[Theta.col_find(tvc)] == tvc: 
     2372                    state = 15 
     2373                else: 
     2374                    state = 14 
     2375 
     2376            elif state == 15: # split out the column v[k] 
     2377                # hh is smallest such that nu[hh] satisfies Lemma 2.25. If it is larger than k+1, 
     2378                # it must be modified, since we are changing that part 
     2379                if k + 1 < hh: 
     2380                    hh = k + 1 
     2381                # hzf is maximal such that indicators line up for nu and zeta 
     2382                if k < hzf: 
     2383                    hzf = k 
     2384                # hb is longest common ancestor of nu and rho 
     2385                if hb >= k: 
     2386                    hb = k 
     2387                    qzb = 0 
     2388                state = 2 
     2389 
     2390            elif state == 16: # backtrack up zeta, updating information about stabilizer vector 
     2391                i = W[k] 
     2392                j = ham_wts[i & 65535] + ham_wts[(i >> 16) & 65535] 
     2393                if j == index and ht == k + 1: ht = k 
     2394                size = size*index 
     2395                index = 0 
     2396                k -= 1 
     2397                state = 13 
     2398 
     2399            elif state == 17: # see if there are any more splits to make from this level of nu (and not zeta) 
     2400                if e[k] == 0: # now is the time to narrow down W[k] by Omega and Phi 
     2401                    # intersect W[k] with each Omega[i] such that v[0]...v[k-1] is in Phi[i] 
     2402                    ii = 0 
     2403                    for i from 0 <= i < k: 
     2404                        ii ^= (1 << v[i]) 
     2405                    for i from 0 <= i <= l: 
     2406                        if Phi[i] & ii == ii: 
     2407                            W[k] &= Omega[i] 
     2408                e[k] = 1 
     2409                # see if there is a vertex to split out 
     2410                i = v[k] + 1 
     2411                while i < ncols and not (1 << i) & W[k]: 
     2412                    i += 1 
     2413                if i < ncols: 
     2414                    v[k] = i 
     2415                    state = 15; continue 
     2416 
     2417                v[k] = -1 
     2418                k -= 1 
     2419                state = 13 
     2420 
     2421            elif state == 18: # the first time nu becomes a discrete partition: set up zeta, our "identity" leaf 
     2422                # initialize counters for zeta: 
     2423                h = k # zeta[h] == nu[h] 
     2424                ht = k # nodes descended from zeta[ht] are all equivalent 
     2425                hzf = k # max such that indicators for zeta and nu agree 
     2426                zeta = PartitionStack(nu) 
     2427                k -= 1 
     2428                rho = PartitionStack(nu) 
     2429                # initialize counters for rho: 
     2430                k_rho = k # number of partitions in rho 
     2431                hzb__h_rho = k # max such that indicators for rho and nu agree - BDM had k+1 
     2432                hb = k # rho[hb] == nu[hb] - BDM had k+1 
     2433                qzb = 0 # Lambda[k] == zb[k], so... 
     2434                state = 13 
     2435 
     2436        # end big while loop 
     2437 
     2438 
     2439 
     2440 

  • sage/graphs/graph_isom.pyx

    r7760 r7762  
    12111211    cdef int *e # 0 or 1, see states 12 and 17 
    12121212    cdef int state # keeps track of place in algorithm 
    1213     cdef int _dig, tvc, tvh, n = G.order() 
     1213    cdef int _dig, tvh, n = G.order() 
    12141214     
    12151215    # trivial case 
    … …  
    16071607            output.append([ Integer(gamma[i]) for i from 0 <= i < n ]) 
    16081608 
    1609             # The variable tvc was set to be the minimum element of W[k] 
     1609            # The variable tvh was set to be the minimum element of W[k] 
    16101610            # the last time we were at state 13 and at a node descending to 
    16111611            # zeta. If this is a minimal cell representative of Theta and 
    … …  
    16141614            # 12, i.e. consider whether we still need to search downward from 
    16151615            # there. TODO: explain why 
    1616             if Theta.elements[tvc] == -1 and lab: ## added "and lab" 
     1616            if Theta.elements[tvh] == -1 and lab: ## added "and lab" 
    16171617                state = 11 
    16181618                continue 
    … …  
    16561656            h = k 
    16571657             
    1658             # set tvc and tvh to the minimum cell representative of W[k] 
     1658            # set tvh to the minimum cell representative of W[k] 
    16591659            # (see states 10 and 14) 
    16601660            for i from 0 <= i < n: 
    16611661                if W[k][i]: 
    1662                     tvc = i 
     1662                    tvh = i 
    16631663                    break 
    1664             tvh = tvc 
    16651664            state = 14 
    16661665             
Note: See TracChangeset for help on using the changeset viewer.