Compute Homology in Python

python homology.py <inputfile>

import capdRedHom
import sys
import ast

def read_dat(filename):
    def read(filename):
        with open(filename, 'r') as input:
            for line in input:
                simplex = map(int, line.split())
                yield simplex

    facets = list(read(filename))
    complex = capdRedHom.SimplicialComplex(facets)
    return complex

def read_pcub(filename):
    def read(filename):
        with open(filename, 'r') as input:
            for line in input:
                cube = ast.literal_eval(line)
                yield cube

    cubes = list(read(filename))
    complex = capdRedHom.CubicalComplex(cubes=cubes)
    return complex

def read(filename):
    if filename.endswith('.dat'):
        return read_dat(filename)
    elif filename.endswith('.pcub_py'):
        return read_pcub(filename)
    else:
        raise RuntimeError("Supported file types: {}".format(['.dat', '.pcub_py']))


if __name__ == '__main__':
    complex = read(sys.argv[1])
    algorithm = capdRedHom.BettiNumbersOverZ(complex)
    betti = algorithm()
    print "Betti: ", betti

Example input files:

cubical torus (.pcub_py extension)

(1, 1),(0, 1),(0, 1),(0, 0)
(0, 0),(0, 1),(0, 0),(0, 1)
(0, 1),(1, 1),(0, 1),(0, 0)
(0, 0),(0, 1),(0, 1),(1, 1)
(0, 1),(0, 0),(1, 1),(0, 1)
(0, 1),(0, 0),(0, 0),(0, 1)
(0, 0),(0, 1),(1, 1),(0, 1)
(1, 1),(0, 1),(1, 1),(0, 1)
(0, 1),(1, 1),(0, 0),(0, 1)
(1, 1),(0, 1),(0, 1),(1, 1)
(0, 1),(1, 1),(1, 1),(0, 1)
(0, 1),(1, 1),(0, 1),(1, 1)
(0, 1),(0, 0),(0, 1),(1, 1)
(1, 1),(0, 1),(0, 0),(0, 1)
(0, 0),(0, 1),(0, 1),(0, 0)
(0, 1),(0, 0),(0, 1),(0, 0)

simplicial torus (.dat extension)

Images Persistent Homology in Python

python image_persistence.py <inputfile>

import sys, os
import gzip
import capdRedHom

try:
    import matplotlib
    matplotlib.use("Agg")
except:
    pass

def read_image(filename):
    if filename.endswith('.gz'):
        fileobj = gzip.open(filename, 'rb')
        filename = filename[:-3]
    else:
        fileobj = open(filename, 'rb')

    if filename.endswith('fits'):
        import pyfits
        data_cube, header_data_cube = pyfits.getdata(fileobj, 0, header=True)
        data = data_cube
    else:
        try:
            from skimage.io import imread
            image = imread(fileobj).astype(float)
            data = image
        except ImportError as ex:
            import PIL.Image
            image = PIL.Image.open(fileobj)
            w, h = image.size
            image_pixels = image.load()
            data = [ [image_pixels[i, j] for j in xrange(w)] for i in xrange(h)]

    return data


def compute_persistence_homology(data):
    algorithm = capdRedHom.persistence.ImagePersistentHomology(data)
    algorithm()
    full_diagram = algorithm.diagram()
    dim_1_diagram = algorithm.diagram(1)
    dim_10_diagram = algorithm.diagram(10) # empty dimension
    return full_diagram, dim_1_diagram, dim_10_diagram


def main(argv=sys.argv):
    image_file = argv[1]
    image_file_bn = os.path.basename(image_file)

    image = read_image(image_file)
    diagram, dim_1_diagram, dim_10_diagram = compute_persistence_homology(image)

    finite, infinite = (len(diagram.finite_intervals), len(diagram.infinite_intervals))
    dim_1_finite, dim_1_infinite = (len(dim_1_diagram.finite_intervals), len(dim_1_diagram.infinite_intervals))
    dim_10_finite, dim_10_infinite = (len(dim_10_diagram.finite_intervals), len(dim_10_diagram.infinite_intervals))

    print "Persistence diagram contains {} finite intervals and {} infinite interals".format(finite, infinite)
    print "Persistence diagram in dimension 1 contains {} finite intervals and {} infinite interals".format(dim_1_finite, dim_1_infinite)
    print "Persistence diagram in dimension 10 contains {} finite intervals and {} infinite interals".format(dim_10_finite, dim_10_infinite)

    with open(image_file_bn + '.intervals', 'w') as out:
        for interval in diagram:
            out.write('{} {}\n'.format(*interval))

    try:
        import matplotlib.pyplot as plt
        plot_diagram = capdRedHom.persistence.PlotDiagram(diagram)
        plot_diagram(plt)
        plt.savefig(image_file_bn + '.png')
        plt.clf()
        matplotlib.pyplot.close("all")
    except ImportError:
        pass


if __name__ == '__main__':
    main(sys.argv)

Example input files:

png

Bottleneck distance in Python

python bottleneck_distance.py <use infinite: 0 1> <diagram1> <diagram2>

import sys
import capdRedHom

def read_diagram(fn):
    with open(fn) as file:
        diagram = capdRedHom.persistence.Diagram.read_from(file)
    return diagram


def main(argv=sys.argv):
    use_inf = bool(argv[1])
    diagram_1 = read_diagram(argv[2])
    diagram_2 = read_diagram(argv[3])

    if capdRedHom.persistence.BottleneckDistance.enabled():
        distance = capdRedHom.persistence.BottleneckDistance(diagram_1, diagram_2, use_inf, 0.001)
        dist = distance()
        print dist
    else:
        print "Not enabled"

if __name__ == '__main__':
    main(sys.argv)

Example input files: diagram1 diagram2

Wasserstein distance in Python

python wasserstein_distance.py <use infinite: 0 1> <diagram1> <diagram2>

import sys
import capdRedHom

def read_diagram(fn):
    with open(fn) as file:
        diagram = capdRedHom.persistence.Diagram.read_from(file)
    return diagram


def main(argv=sys.argv):
    use_inf = bool(argv[1])
    diagram_1 = read_diagram(argv[2])
    diagram_2 = read_diagram(argv[3])

    if capdRedHom.persistence.WassersteinDistance.enabled():
        distance = capdRedHom.persistence.WassersteinDistance(diagram_1, diagram_2, use_inf)
        dist = distance()
        print dist
    else:
        print "Not enabled"

if __name__ == '__main__':
    main(sys.argv)

Example input files: diagram1 diagram2