In this notebook we will use Polymake to compute the persistent homology of a Vietoris-Rips complex.

A popular book about these topics is Computational Topology: An Introduction by Herbert Edelsbrunner and John Harer. A shorter introduction with a focus on the algorithms involved can be found in the paper Computing Persistent Homology by Afra Zomorodian and Gunnar Carlsson.

After computing persistent homology we will visualize it as a barcode diagram using Julias "Plots" package.

In [1]:
versioninfo()
Julia Version 1.4.0
Commit b8e9a9ecc6 (2020-03-21 16:36 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, ivybridge)
In [1]:
import Pkg
Pkg.activate(".")
Pkg.status()

using Polymake
using Plots 
using Distances
 Activating environment at `~/Desktop/PMJLFork/Polymake.jl/examples/persistent_homology/Project.toml`
Status `~/Desktop/PMJLFork/Polymake.jl/examples/persistent_homology/Project.toml`
  [1f15a43c] CxxWrap v0.10.1
  [b4f34e82] Distances v0.9.0
  [91a5bcdd] Plots v1.3.6
  [d720cf60] Polymake v0.4.2 [`~/.julia/dev/Polymake`]
polymake version 4.0
Copyright (c) 1997-2020
Ewgenij Gawrilow, Michael Joswig, and the polymake team
Technische Universität Berlin, Germany
https://polymake.org

This is free software licensed under GPL; see the source for copying conditions.
There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

We wil begin by randomly generating a set of points and plotting it, if it is in dimension two or three.

In [2]:
n = 8
dim = 3

points = rand(Float64, (n,dim))

fig = plot()

if dim == 2
    plot!(points[:,1],points[:,2],
    seriestype=:scatter, markersize = 5,legend=false,aspect_ratio =:equal)
    @show fig
end

if dim == 3
    plot!(points[:,1],points[:,2], points[:,3],
        seriestype=:scatter, markersize = 5,legend=false,aspect_ratio =:equal)
    @show fig
end
fig = Plot{Plots.GRBackend() n=1}
Out[2]:

Let $V = \{v_1,...,v_n\}$ be our vertex set. The Vietoris Rips complex of $V$ with respect to some $\epsilon > 0$ is defined as

$\operatorname{VR}_\epsilon(V) := \{\sigma \subset V \mid \operatorname{dist}(v_i,v_j)\leq \epsilon$, $v_i,v_j \in \sigma\}$

Polymake requires a distance matrix to compute the Vietoris rips complex. Here we use the "Distances" package to compute a $n x n$ matrix of pairwise distances, where $n$ is the number of points. We can specify which metric we want to use and choose the Euclidean metric.

Setting dims = 1 specifies that we want the pairwise distances of the rows of our matrix points.

In [3]:
distances = Distances.pairwise(Distances.Euclidean(),points,dims=1)
Out[3]:
8×8 Array{Float64,2}:
 0.0       0.787939  0.461256  0.949579  …  0.892467  0.430361  0.477265
 0.787939  0.0       0.369462  0.296285     0.625871  0.797832  0.757861
 0.461256  0.369462  0.0       0.551528     0.552891  0.632143  0.482414
 0.949579  0.296285  0.551528  0.0          0.546088  0.895457  1.00846
 0.341375  0.545702  0.316362  0.782494     0.86664   0.512323  0.35958
 0.892467  0.625871  0.552891  0.546088  …  0.0       1.03158   0.940508
 0.430361  0.797832  0.632143  0.895457     1.03158   0.0       0.823685
 0.477265  0.757861  0.482414  1.00846      0.940508  0.823685  0.0

To get a filtration we construct Vietoris-Rips complexes with values $\epsilon_i$ which increase by a specified stepwidth in each step. These steps are called frames.

Since all points in our vertex set appear for every $\epsilon > 0$ we need to decide in which frames they should appear in a different manner. To this end we define a Polymake.Array() which stores the frame of appearance for each point in vertex set $V$. A natural choice is for all to appear in the zeroth frame or in consecutive frames starting at $0$.

In [4]:
a = Polymake.Array(zeros(Int8,n))
Out[4]:
pm::Array<pm::Integer>
0 0 0 0 0 0 0 0

Now we can build the actual filtration via Polymakes topaz module. Along our distance matrix and the specified degrees for the points we pass a stepwidth and another Integer. This Integer specifies the dimension up to which the VR-complex is computed in each step.

Afterwards we use anoter function of the topaz module to compute the persistent homology of our filtered complex.

The output is an array of a list of pairs. Each list represents the homology classes of one dimension in ascending order. Each pair represents one homology class. The values stored in the pairs specify the frame in which the class is born and the frame in which it dies. A value of $-1$ indicates that the homology class is not killed at any point.

In [5]:
filtration = @pm topaz.vietoris_rips_filtration{Rational}(distances, a, 0.1, 4)
pershom = topaz.persistent_homology(filtration)
Out[5]:
pm::Array<std::list<std::pair<long, long>, std::allocator<std::pair<long, long> > >>
{(0 3) (0 4) (0 4) (0 4) (0 4) (0 5) (0 6) (0 -1)}
{}
{}
{}
{(9 -1) (10 -1) (10 -1) (10 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1) (11 -1)}

Now we extract the highest frame number in which a homology class dies. We need this for our visualization later on.

In [6]:
inf = 0
for i in 1:length(pershom)
    if(length(pershom[i])>0)
        inf = max(inf,maximum(maximum(pershom[i])))
    end
end
inf += 2
Out[6]:
13

The next code block converts the information of our pershom object into different Julia arrays we will need to plot the diagram. We need different containers for:

  • the points we want to plot, inidcating birth and death of homology classes: x and y
  • lines between points, where each line represents the lifetime of a homology class: lines
  • arrows from a point of birht to infinity: arrows
  • some purely informative stuff: breaklines, custom_yticks and custom_ylabels
In [7]:
x = []
y = []
lines = []
arrows = []
levels = [0]
breaklines = []
custom_yticks = [0.0]
custom_ylabels = []

level = 1
for (index,dims) in enumerate(pershom)
    for elements in dims
        append!(x,[elements[1]])
        append!(y,[level])
        
        if(elements[2] == -1)
            append!(arrows,[[[elements[1], inf], [level, level]]])
        else 
            append!(x,[elements[2]])
            append!(y,[level])
            append!(lines,[[[elements[1], elements[2]], [level, level]]])
        end
        level+=1
    end
    
    append!(levels,level)
    append!(breaklines,[[[-0.25,inf],[level, level]]])
    append!(custom_ylabels, ["H$(index-1)"])
    level+=1
end

for i in 2:length(levels)
    append!(custom_yticks, levels[i-1] + (levels[i]-levels[i-1])/2)
end

Now we can use the generated information to generate the bardiagram corresponding to our computations. This concludes our tutorial.

In [8]:
fig = plot()
for line in lines    
    plot!(line[1], line[2], linecolor=:steelblue)
end
for arrow in arrows
    plot!(arrow[1],arrow[2], arrow = true, linecolor=:steelblue)
end
for breakline in breaklines[1:size(breaklines,1)-1]
    plot!(breakline[1],breakline[2], line =:dash, linecolor=:black)
end
scatter!(x,y,xlims = (-0.25,inf),xticks = 0:1:inf-1, yticks = (custom_yticks[2:size(custom_yticks,1)],custom_ylabels), xlabel="Frame",legend=false)
fig
Out[8]:

Another popular way of visualizing persistence is to draw a diagram where the $x$ and $y$-axis represent frames and we draw a point $(\operatorname{frame of birth}, \operatorname{frame of death})$ for each homology clas. The homology classes correpsonding to different dimensions are colored differently.

In [35]:
x = []
y = []
fig = plot()
plot!([-0.25,inf],[-0.25,inf], linecolor=:black,label = "")
plot!([-0.25,inf],[inf,inf], line=:dash, linecolor=:black, label = "")

for (id,dims) in enumerate(pershom)
    for (jd,elements) in enumerate(dims)
        #Splitting the scatter draws to only get legend of initial point in each dimension. Very hacky.
        if jd == 1
            if elements[2] == -1
                scatter!([elements.first],[inf],xlims = (-0.25,inf),ylims = (0,inf+1),xticks = 0:1:inf-1, 
                yticks = 0:1:inf-1,label = string(id-1),legend =:bottomright, color=id)
            else
                scatter!([elements.first],[elements.second],xlims = (-0.25,inf),ylims = (0,inf+1),xticks = 0:1:inf-1, 
                yticks = 0:1:inf-1,label = string(id-1),legend =:bottomright, color=id)
            end
        else
            append!(x,elements.first)
            if elements[2] == -1
                append!(y,Int(inf))
            else
                append!(y,Int(elements.second))
            end
        end
        
        
    scatter!(x,y,xlims = (-0.25,inf),ylims = (0,inf+1),xticks = 0:1:inf-1, 
           yticks = 0:1:inf-1,label = "",  color=id)
    
    x = []
    y = []
    end
end

fig
Out[35]:
In [ ]: