Computing the Homology of a Simplicial Complex

homology.py

import numpy as np
import functools as func
import itertools as itertools

# spcx::Set<Tuple> --> int
def vertex_count(spcx):
  # Reduce Set to max value amongst Tuples in the set
  return func.reduce(lambda acc, e: max(acc, max(e)), spcx, 0) + 1

# spcx::Set<Tuple>, n::int --> Dict<Tuple, int>
def C_n_basis(scpx, n):
  # Construct 
  basis = dict()
  i = 0
  for face in scpx:
    for f in itertools.combinations(face, n + 1):
      if (f not in basis):
        basis[f] = i
        i += 1
  return basis

# spcx::Set<Tuple>, n::int --> Set<Tuple>
def C_n(spcx, n):
  basis = set()
  for face in spcx:
    for f in itertools.combinations(face, n + 1):
      # Set union with reassignment
      basis |= {f}
  return basis

# Given a face, return its boundary in respect to the given basis
# face::Tuple, basis::Dict<Tuple, int> --> int[]
def boundary_vec(face, basis):
  vec = [0] * len(basis)
  for i in range(0, len(face)):
    new_vec = list(face)
    del new_vec[i] # remove i-th element of face
    index = basis[tuple(new_vec)]
    vec[index] = (-1) ** i
  return vec

# scpx::Set<Tuple>, n::int --> np.matrix
# assumes n > 0
def boundary_matrix(scpx, n):
  colBasis = C_n(scpx, n)
  rowBasis = C_n_basis(scpx, n - 1)
  cols = []
  for c in colBasis:
    cols.append(boundary_vec(c, rowBasis))
  # Construct numpy matrix column by column
  return np.column_stack(cols)

# Calculate the dimension of n-th homology
# scpx::Set<Tuple>, n::int --> int
def homology_n(scpx, n):
  dimIm_d_n1 = 0
  dimKer_d_n = 0
  # Don't bother calculating dimIm_d_3 because we know it's 0
  if (n != 2):
    d_n1 = boundary_matrix(C_n(scpx, n + 1), n + 1)
    dimIm_d_n1 = np.linalg.matrix_rank(d_n1)
  if (n == 0):
      # Don't bother calculating boundary matrix for the case n = 0
      # It's just the number of vertices
    dimKer_d_n = vertex_count(scpx)
  else:
    d_n = boundary_matrix(C_n(scpx, n), n)
    # Rank-Nullity Theorem
    dimKer_d_n = d_n.shape[1] - np.linalg.matrix_rank(d_n)
  # Calculate dimension of n-th homology
  return dimKer_d_n - dimIm_d_n1

# Print the Homology of a Simplicial Complex
def homology(name, scpx):
  print(f'{name}')
  print(f'{"-"*len(name)}')
  print(f'dim H_0: {homology_n(scpx, 0)}')
  print(f'dim H_1: {homology_n(scpx, 1)}')
  print(f'dim H_2: {homology_n(scpx, 2)}')
  print('\n')

# Top level simplicies
# Cannot be a Set of Sets because Python Sets are non-hashable types
#  and elements of a Python Set must be hashable. Tuples are hashable
# Tuples must be in increasing order

sphere = {
  (0,1,5), (0,3,5), (2,3,5), (1,2,5),
  (0,1,4), (0,3,4), (2,3,4), (1,2,4)
}

sphere_torus_d = {
  (0,2,7), (0,3,7), (3,7,8), (3,4,8), (2,4,8), (0,2,4),
  (1,2,5), (2,5,7), (5,6,7), (6,7,8), (1,6,8), (1,2,8), 
  (0,1,3), (1,3,5), (3,4,5), (4,5,6), (0,4,6), (0,1,6),
  # sphere points now --- disconnected
  (9,10,14), (9,12,14), (11,12,14), (10,11,14),
  (9,10,13), (9,12,13), (11,12,13), (10,11,13)
}

sphere_torus_c = {
  (0,2,7), (0,3,7), (3,7,8), (3,4,8), (2,4,8), (0,2,4),
  (1,2,5), (2,5,7), (5,6,7), (6,7,8), (1,6,8), (1,2,8), 
  (0,1,3), (1,3,5), (3,4,5), (4,5,6), (0,4,6), (0,1,6),
  # sphere points now --- connected
  (8,9,13), (8,11,13), (10,11,13), (9,10,13),
  (8,9,12), (8,11,12), (10,11,12), (9,10,12)
}

torus = {
  (0,2,7), (0,3,7), (3,7,8), (3,4,8), (2,4,8), (0,2,4),
  (1,2,5), (2,5,7), (5,6,7), (6,7,8), (1,6,8), (1,2,8), 
  (0,1,3), (1,3,5), (3,4,5), (4,5,6), (0,4,6), (0,1,6)
}

real_proj = {
  (0,1,4), (1,4,5), (1,2,5),
  (0,2,5), (0,3,5), (0,1,3), (1,2,3),
  (2,3,4), (0,2,4), (3,4,5)
}

homology("Sphere+Torus Disconnected", sphere_torus_d)
homology("Sphere+Torus Connected", sphere_torus_c)
homology("Sphere", sphere)
homology("Torus", torus)
homology("RP^2", real_proj)