Finding nested isotropic spaces in characteristic 2.

isotropic.py

from itertools import combinations, product
import numpy as np
import csv

def get_lucas_matrix(genus=3):
    # generate V, all non-zero binary 2*genus tuples.
    V = (np.array(m) 
         for m in product(*( [range(2)]*genus*2 ))
         if np.array(m).sum()>0)

    # generate all distinct unordered g-tuples of V
    genuslets = combinations(V, genus)
    isotropic_tuples = (t for t in genuslets
                        if is_isotropic(t))

    smalls = set()
    bigs = set()
    
    for vectors in isotropic_tuples:
        space = Space(np.array(vectors))
        if space.rank == genus:
            bigs.add(space)
        elif space.rank == genus-1:
            smalls.add(space)

    return np.array([[small < big for big in bigs] 
                     for small in smalls], 
                    dtype=int)

def has_rank_2(a,b,c):
    return ((a+b+c)%2).sum() == 0

def is_isotropic(g_tuple):
    return all(map(orthogonal, 
                   combinations(g_tuple,2)))

# # defines the bilinear form
# OMEGA = np.array([[0,1,0,0,0,0],
#                   [1,0,0,0,0,0],
#                   [0,0,0,1,0,0],
#                   [0,0,1,0,0,0],
#                   [0,0,0,0,0,1],
#                   [0,0,0,0,1,0]])

# # def orthogonal((x,y)):
# #     return x.T.dot(OMEGA).dot(y) % 2 == 0

def orthogonal((x,y)):
    return x[::-1].dot(y) %2 == 0

def span(a,b,c):
    return set(map(tuple,
                   [a,b,c,
                    (a+b)%2, (a+c)%2, (b+c)%2,
                    (a+b+c)%2]))

def mod2rank(M, in_place=False):
  M = np.array(M, copy=not in_place)
  for index, column in enumerate(M.T):
    nonzero_rows = np.flatnonzero(column)
    if any(nonzero_rows):
      first_one, other_ones = ( nonzero_rows[0], 
                                nonzero_rows[1:] )
      M[other_ones] = (M[other_ones]+M[first_one])%2
      M[first_one, index+1:] = 0
  return M.sum()

class Space:
    def __init__(self,M):
        M = np.array(M) #, copy=False)
        pivots = []
        rank = 0

        # for index,column in enumerate(M.T):
        #     nonzeros = np.flatnonzero(column)
        #     if len(nonzeros)>1:
        #         a, others = nonzeros[0], nonzeros[1:]
        #         M[others] = (M[a] + M[others])%2
        #         try:
        #             M[a], M[rank] = M[rank], M[a]
        #         except:
        #             import pdb; pdb.set_trace()
        #     if len(nonzeros)>0:
        #         pivots.append(index)
        #         rank += 1

        for col in range(len(M[0])):
            for row in range(rank,len(M)):
                if M[row][col] == 1:
                    for s in range(len(M)):
                        if M[s][col]==1 and (s!=row):
                            M[s] = (M[row] + M[s]) %2
                    M[row],M[rank]=M[rank],M[row]
                    pivots.append(col)
                    rank += 1
                    break
        self.rank = rank                    
        self.basis = tuple(map(tuple, M[:self.rank]))
        self.lastpivot = pivots[-1]

    def dim(self):
        return len(self.basis)

    def __gt__(self, N):
        # sorry for this assumption
        assert N.rank + 1 == self.rank

        for vector in self.basis:
            if Space(list(N.basis)+[vector]) == self:
                return True
        return False

    def __lt__(self, N):
        return N > self

    def __eq__(self, N):
        return self.__hash__() == N.__hash__()

    def __hash__(self):
        return self.basis.__hash__()

if __name__=="__main__":
    data = get_lucas_matrix()
    # with open("lucas_matrix.csv",'w') as f:
    #     w = csv.writer(f)
    #     for d in data:
    #         w.writerow(d)

    # check each 2D space is in exactly 3 3D spaces.
    assert all(data.sum(1) == 3)
    print mod2rank(data, in_place=True)

This is currently broken... row-reduction should be made pythonic.