Computing_Homology.py

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Sun Dec  2 15:29:09 2018

@author: Jeremy Kun
https://github.com/j2kun/computing-homology/blob/master/homology.py
"""

import numpy
import numpy.linalg

def rowSwap(A, i, j):
   temp = numpy.copy(A[i, :])
   A[i, :] = A[j, :]
   A[j, :] = temp

def colSwap(A, i, j):
   temp = numpy.copy(A[:, i])
   A[:, i] = A[:, j]
   A[:, j] = temp

def scaleCol(A, i, c):
   A[:, i] *= c*numpy.ones(A.shape[0])

def scaleRow(A, i, c):
   A[i, :] *= c*numpy.ones(A.shape[1])

def colCombine(A, addTo, scaleCol, scaleAmt):
   A[:, addTo] += scaleAmt * A[:, scaleCol]

def rowCombine(A, addTo, scaleRow, scaleAmt):
   A[addTo, :] += scaleAmt * A[scaleRow, :]

def simultaneousReduce(A, B):
   if A.shape[1] != B.shape[0]:
      raise Exception("Matrices have the wrong shape.")

   numRows, numCols = A.shape

   i,j = 0,0
   while True:
      if i >= numRows or j >= numCols:
         break

      if A[i,j] == 0:
         nonzeroCol = j
         while nonzeroCol < numCols and A[i,nonzeroCol] == 0:
            nonzeroCol += 1

         if nonzeroCol == numCols:
            i += 1
            continue

         colSwap(A, j, nonzeroCol)
         rowSwap(B, j, nonzeroCol)

      pivot = A[i,j]
      scaleCol(A, j, 1.0 / pivot)
      scaleRow(B, j, 1.0 / pivot)

      for otherCol in range(0, numCols):
         if otherCol == j:
            continue
         if A[i, otherCol] != 0:
            scaleAmt = -A[i, otherCol]
            colCombine(A, otherCol, j, scaleAmt)
            rowCombine(B, j, otherCol, -scaleAmt)

      i += 1; j+= 1

   return A,B


def finishRowReducing(B):
   numRows, numCols = B.shape

   i,j = 0,0
   while True:
      if i >= numRows or j >= numCols:
         break

      if B[i, j] == 0:
         nonzeroRow = i
         while nonzeroRow < numRows and B[nonzeroRow, j] == 0:
            nonzeroRow += 1

         if nonzeroRow == numRows:
            j += 1
            continue

         rowSwap(B, i, nonzeroRow)

      pivot = B[i, j]
      scaleRow(B, i, 1.0 / pivot)

      for otherRow in range(0, numRows):
         if otherRow == i:
            continue
         if B[otherRow, j] != 0:
            scaleAmt = -B[otherRow, j]
            rowCombine(B, otherRow, i, scaleAmt)

      i += 1; j+= 1

   return B


def numPivotCols(A):
   z = numpy.zeros(A.shape[0])
   return [numpy.all(A[:, j] == z) for j in range(A.shape[1])].count(False)


def numPivotRows(A):
   z = numpy.zeros(A.shape[1])
   return [numpy.all(A[i, :] == z) for i in range(A.shape[0])].count(False)


def bettiNumber(d_k, d_kplus1):
   A, B = numpy.copy(d_k), numpy.copy(d_kplus1)
   simultaneousReduce(A, B)
   finishRowReducing(B)

   dimKChains = A.shape[1]
   print(dimKChains)
   kernelDim = dimKChains - numPivotCols(A)
   print(kernelDim)
   imageDim = numPivotRows(B)
   print(imageDim)

   return kernelDim - imageDim


bd0 = numpy.array([[0.,0,0,0,0]])
bd1 = numpy.array([[-1.,-1,-1,-1,0,0,0,0], [1,0,0,0,-1,-1,0,0],
     [0,1,0,0,1,0,-1,-1], [0,0,1,0,0,1,1,0], [0,0,0,1,0,0,0,1]])
bd2 = numpy.array([[1.,1,0,0],[-1,0,1,0],[0,-1,-1,0],
     [0,0,0,0],[1,0,0,1],[0,1,0,-1],
     [0,0,1,1],[0,0,0,0]])
bd3 = numpy.array([[-1.],[1],[-1],[1]])


print("Example complex from post")
print("0th homology: %d" % bettiNumber(bd0,bd1))
print("1st homology: %d" % bettiNumber(bd1,bd2))
print("2nd homology: %d" % bettiNumber(bd2,bd3))


mobiusD1 = numpy.array([
   [-1.,-1,-1,-1, 0, 0, 0, 0, 0, 0],
   [ 1, 0, 0, 0,-1,-1,-1, 0, 0, 0],
   [ 0, 1, 0, 0, 1, 0, 0,-1,-1, 0],
   [ 0, 0, 0, 1, 0, 0, 1, 0, 1, 1],
])

mobiusD2 = numpy.array([
   [ 1., 0, 0, 0, 1],
   [ 0, 0, 0, 1, 0],
   [-1, 0, 0, 0, 0],
   [ 0, 0, 0,-1,-1],
   [ 0, 1, 0, 0, 0],
   [ 1,-1, 0, 0, 0],
   [ 0, 0, 0, 0, 1],
   [ 0, 1, 1, 0, 0],
   [ 0, 0,-1, 1, 0],
   [ 0, 0, 1, 0, 0],
])

print("Mobius Band")
print("1st homology: %d" % bettiNumber(mobiusD1, mobiusD2))