{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook we will use [Polymake](https://polymake.org/doku.php) to compute the persistent homology of a Vietoris-Rips complex. \n",
"\n",
"A popular book about these topics is [Computational Topology: An Introduction](https://www.researchgate.net/publication/220692408_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](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.10.5064) by Afra Zomorodian and Gunnar Carlsson. \n",
"\n",
"After computing persistent homology we will visualize it as a barcode diagram using [Julias](https://julialang.org/) \"Plots\" package."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Julia Version 1.4.0\n",
"Commit b8e9a9ecc6 (2020-03-21 16:36 UTC)\n",
"Platform Info:\n",
" OS: Linux (x86_64-pc-linux-gnu)\n",
" CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz\n",
" WORD_SIZE: 64\n",
" LIBM: libopenlibm\n",
" LLVM: libLLVM-8.0.1 (ORCJIT, ivybridge)\n"
]
}
],
"source": [
"versioninfo()"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m environment at `~/Desktop/PMJLFork/Polymake.jl/examples/persistent_homology/Project.toml`\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1mStatus\u001b[22m\u001b[39m `~/Desktop/PMJLFork/Polymake.jl/examples/persistent_homology/Project.toml`\n",
" \u001b[90m [1f15a43c]\u001b[39m\u001b[37m CxxWrap v0.10.1\u001b[39m\n",
" \u001b[90m [b4f34e82]\u001b[39m\u001b[37m Distances v0.9.0\u001b[39m\n",
" \u001b[90m [91a5bcdd]\u001b[39m\u001b[37m Plots v1.3.6\u001b[39m\n",
" \u001b[90m [d720cf60]\u001b[39m\u001b[37m Polymake v0.4.2 [`~/.julia/dev/Polymake`]\u001b[39m\n",
"polymake version 4.0\n",
"Copyright (c) 1997-2020\n",
"Ewgenij Gawrilow, Michael Joswig, and the polymake team\n",
"Technische Universität Berlin, Germany\n",
"https://polymake.org\n",
"\n",
"This is free software licensed under GPL; see the source for copying conditions.\n",
"There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\n",
"\n"
]
}
],
"source": [
"import Pkg\n",
"Pkg.activate(\".\")\n",
"Pkg.status()\n",
"\n",
"using Polymake\n",
"using Plots \n",
"using Distances"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We wil begin by randomly generating a set of points and plotting it, if it is in dimension two or three."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"fig = Plot{Plots.GRBackend() n=1}\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 8\n",
"dim = 3\n",
"\n",
"points = rand(Float64, (n,dim))\n",
"\n",
"fig = plot()\n",
"\n",
"if dim == 2\n",
" plot!(points[:,1],points[:,2],\n",
" seriestype=:scatter, markersize = 5,legend=false,aspect_ratio =:equal)\n",
" @show fig\n",
"end\n",
"\n",
"if dim == 3\n",
" plot!(points[:,1],points[:,2], points[:,3],\n",
" seriestype=:scatter, markersize = 5,legend=false,aspect_ratio =:equal)\n",
" @show fig\n",
"end\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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\n",
"\n",
"$\\operatorname{VR}_\\epsilon(V) := \\{\\sigma \\subset V \\mid \\operatorname{dist}(v_i,v_j)\\leq \\epsilon$, $v_i,v_j \\in \\sigma\\}$\n",
"\n",
"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. \n",
"We can specify which metric we want to use and choose the Euclidean metric. \n",
"\n",
"Setting `dims = 1` specifies that we want the pairwise distances of the rows of our matrix `points`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"8×8 Array{Float64,2}:\n",
" 0.0 0.787939 0.461256 0.949579 … 0.892467 0.430361 0.477265\n",
" 0.787939 0.0 0.369462 0.296285 0.625871 0.797832 0.757861\n",
" 0.461256 0.369462 0.0 0.551528 0.552891 0.632143 0.482414\n",
" 0.949579 0.296285 0.551528 0.0 0.546088 0.895457 1.00846\n",
" 0.341375 0.545702 0.316362 0.782494 0.86664 0.512323 0.35958\n",
" 0.892467 0.625871 0.552891 0.546088 … 0.0 1.03158 0.940508\n",
" 0.430361 0.797832 0.632143 0.895457 1.03158 0.0 0.823685\n",
" 0.477265 0.757861 0.482414 1.00846 0.940508 0.823685 0.0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"distances = Distances.pairwise(Distances.Euclidean(),points,dims=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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__. \n",
"\n",
"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$. \n",
"A natural choice is for all to appear in the zeroth frame or in consecutive frames starting at $0$. \n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pm::Array\n",
"0 0 0 0 0 0 0 0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = Polymake.Array(zeros(Int8,n))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"This Integer specifies the dimension up to which the VR-complex is computed in each step. \n",
"\n",
"Afterwards we use anoter function of the __topaz__ module to compute the persistent homology of our filtered complex. \n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"pm::Array, std::allocator > >>\n",
"{(0 3) (0 4) (0 4) (0 4) (0 4) (0 5) (0 6) (0 -1)}\n",
"{}\n",
"{}\n",
"{}\n",
"{(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)}\n"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"filtration = @pm topaz.vietoris_rips_filtration{Rational}(distances, a, 0.1, 4)\n",
"pershom = topaz.persistent_homology(filtration)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we extract the highest frame number in which a homology class dies. We need this for our visualization later on."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"13"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inf = 0\n",
"for i in 1:length(pershom)\n",
" if(length(pershom[i])>0)\n",
" inf = max(inf,maximum(maximum(pershom[i])))\n",
" end\n",
"end\n",
"inf += 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:\n",
"* the points we want to plot, inidcating birth and death of homology classes: `x` and `y`\n",
"* lines between points, where each line represents the lifetime of a homology class: `lines`\n",
"* arrows from a point of birht to infinity: `arrows`\n",
"* some purely informative stuff: `breaklines`, `custom_yticks` and `custom_ylabels`"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"x = []\n",
"y = []\n",
"lines = []\n",
"arrows = []\n",
"levels = [0]\n",
"breaklines = []\n",
"custom_yticks = [0.0]\n",
"custom_ylabels = []\n",
"\n",
"level = 1\n",
"for (index,dims) in enumerate(pershom)\n",
" for elements in dims\n",
" append!(x,[elements[1]])\n",
" append!(y,[level])\n",
" \n",
" if(elements[2] == -1)\n",
" append!(arrows,[[[elements[1], inf], [level, level]]])\n",
" else \n",
" append!(x,[elements[2]])\n",
" append!(y,[level])\n",
" append!(lines,[[[elements[1], elements[2]], [level, level]]])\n",
" end\n",
" level+=1\n",
" end\n",
" \n",
" append!(levels,level)\n",
" append!(breaklines,[[[-0.25,inf],[level, level]]])\n",
" append!(custom_ylabels, [\"H$(index-1)\"])\n",
" level+=1\n",
"end\n",
"\n",
"for i in 2:length(levels)\n",
" append!(custom_yticks, levels[i-1] + (levels[i]-levels[i-1])/2)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use the generated information to generate the bardiagram corresponding to our computations. This concludes our tutorial."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fig = plot()\n",
"for line in lines \n",
" plot!(line[1], line[2], linecolor=:steelblue)\n",
"end\n",
"for arrow in arrows\n",
" plot!(arrow[1],arrow[2], arrow = true, linecolor=:steelblue)\n",
"end\n",
"for breakline in breaklines[1:size(breaklines,1)-1]\n",
" plot!(breakline[1],breakline[2], line =:dash, linecolor=:black)\n",
"end\n",
"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)\n",
"fig"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = []\n",
"y = []\n",
"fig = plot()\n",
"plot!([-0.25,inf],[-0.25,inf], linecolor=:black,label = \"\")\n",
"plot!([-0.25,inf],[inf,inf], line=:dash, linecolor=:black, label = \"\")\n",
"\n",
"for (id,dims) in enumerate(pershom)\n",
" for (jd,elements) in enumerate(dims)\n",
" #Splitting the scatter draws to only get legend of initial point in each dimension. Very hacky.\n",
" if jd == 1\n",
" if elements[2] == -1\n",
" scatter!([elements.first],[inf],xlims = (-0.25,inf),ylims = (0,inf+1),xticks = 0:1:inf-1, \n",
" yticks = 0:1:inf-1,label = string(id-1),legend =:bottomright, color=id)\n",
" else\n",
" scatter!([elements.first],[elements.second],xlims = (-0.25,inf),ylims = (0,inf+1),xticks = 0:1:inf-1, \n",
" yticks = 0:1:inf-1,label = string(id-1),legend =:bottomright, color=id)\n",
" end\n",
" else\n",
" append!(x,elements.first)\n",
" if elements[2] == -1\n",
" append!(y,Int(inf))\n",
" else\n",
" append!(y,Int(elements.second))\n",
" end\n",
" end\n",
" \n",
" \n",
" scatter!(x,y,xlims = (-0.25,inf),ylims = (0,inf+1),xticks = 0:1:inf-1, \n",
" yticks = 0:1:inf-1,label = \"\", color=id)\n",
" \n",
" x = []\n",
" y = []\n",
" end\n",
"end\n",
"\n",
"fig"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.4.0",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}