Compute the Betti numbers of a graph over Z/2Z

betti_bin.py

# Compute the Betti numbers of a graph over Z/2Z.
# 
# If G is a networkx graph then set
#   C = networkx.find_cliques(G)
# This enumerates maximal cliques.
#
# Pass C to the function betti_bin,
# which returns a list of Betti numbers.
#
# RA, 2017-11-08 (CC-BY-4.0)
#
# Ref: 
#   A. Zomorodian, Computational topology (Notes), 2009
#   http://www.ams.org/meetings/short-courses/zomorodian-notes.pdf
#
# The computation of the rank is adapted from
#   https://triangleinequality.wordpress.com/2014/01/23/computing-homology/
# see also
#   https://gist.github.com/numpde/9584779ad235c6ee19be7a6bb87e8af5
#
# Gist for this function:
#   https://gist.github.com/numpde/a3f3476e601041da3cb968f31e95409d
#
def betti_bin(C) :
    from itertools import combinations as subcliques
    
    # Each maximal clique as a sorted tuple
    C = [tuple(sorted(mc)) for mc in C]
    
    # Betti numbers
    B = []
    
    # (k-1)-chain group
    Sl = dict()
    
    # Iterate over the dimension
    for k in range(0, max(len(s) for s in C)) :
        
        # Get all (k+1)-cliques, i.e. k-simplices, from max cliques mc
        Sk = set(c for mc in C for c in subcliques(mc, k+1))
        # Check that each simplex is in increasing order
        assert(all((list(s) == sorted(s)) for s in Sk))
        # Assign an ID to each simplex, in lexicographic order
        # (This ordering makes subsequent computations faster)
        Sk = dict(zip(sorted(Sk), range(0, len(Sk))))
        
        # Sk is now a representation of the k-chain group
        
        # "Default" Betti number (if rank is 0)
        B.append(len(Sk))
        
        # J2I is a mapped representation of the boundary operator
        # (Understood to have boolean entries instead of +1 / -1)
        
        J2I = {
            j : set(Sl[s] for s in subcliques(ks, k) if s)
            for (ks, j) in Sk.items()
        }
        
        Sl = Sk
        
        # Transpose J2I
        
        I2J = {
            i : set()
            for I in J2I.values() for i in I
        }
        
        for (j, I) in J2I.items() :
            for i in I :
                I2J[i].add(j)
                
        # Compute the rank of the boundary operator
        # The rank = The # of completed while loops
        # It's usually the most time-consuming part
        while J2I and I2J :
            # Default pivot option
            I = next(iter(J2I.values()))
            J = I2J[next(iter(I))]
            
            # An alternative option
            JA = next(iter(I2J.values()))
            IA = J2I[next(iter(JA))]
            
            # Choose the option with fewer indices
            # (reducing the number of loops below)
            if ((len(IA) + len(JA)) < (len(I) + len(J))) : (I, J) = (IA, JA)

            for i in I : 
                I2J[i] = J.symmetric_difference(I2J[i])
                if not I2J[i] : del I2J[i]

            for j in J : 
                J2I[j] = I.symmetric_difference(J2I[j])
                if not J2I[j] : del J2I[j]
            
            # Let
            #   D_k : C_k --> C_{k-1}
            # be the boundary operator. Recall that
            #   H_k = ker D_k / im D_{k+1}

            # Rank of  D_k  increases, therefore:
            B[k-1] -= 1
            # Nullspace dimension of  D_k  decreases, therefore:
            B[k]   -= 1
    
    # Remove trailing zeros
    while B and (B[-1] == 0) : B.pop()
    
    return B