{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook we will demonstrate how Polymake can be used to compute a Vietoris Rips complex and its homology on a random two dimensional point cloud. To read up on these concepts we suggest to take a look at: [Computational Topology: An Introduction](https://www.researchgate.net/publication/220692408_Computational_Topology_An_Introduction) by Herbert Edelsbrunner and John Harer. \n",
"\n",
"We will begin by checking on our system by calling `versioninfo()`. We can see the version of our Julia Kernel and some system specs."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"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": "markdown",
"metadata": {},
"source": [
"Then we check the status of our package manager and tell Julia which packages to use."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[32m\u001b[1mStatus\u001b[22m\u001b[39m `~/.julia/environments/v1.4/Project.toml`\n",
" \u001b[90m [1f15a43c]\u001b[39m\u001b[37m CxxWrap v0.9.1\u001b[39m\n",
" \u001b[90m [b4f34e82]\u001b[39m\u001b[37m Distances v0.8.2\u001b[39m\n",
" \u001b[90m [7073ff75]\u001b[39m\u001b[37m IJulia v1.21.1\u001b[39m\n",
" \u001b[90m [91a5bcdd]\u001b[39m\u001b[37m Plots v1.0.5\u001b[39m\n",
" \u001b[90m [d720cf60]\u001b[39m\u001b[37m Polymake v0.4.1\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": [
"using Pkg\n",
"Pkg.status()\n",
"\n",
"using Polymake\n",
"using Plots \n",
"using Distances"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With this our setup is complete. We will begin by generating a random set of points via some of Julia's basic functionality. Then we will plot these points using the \"Plots\" package."
]
},
{
"cell_type": "code",
"execution_count": 117,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"16×2 Array{Float64,2}:\n",
" 0.073776 0.861983\n",
" 0.255749 0.216406\n",
" 0.837399 0.18367\n",
" 0.451908 0.482145\n",
" 0.968421 0.898811\n",
" 0.727573 0.698121\n",
" 0.176284 0.384034\n",
" 0.0737789 0.301895\n",
" 0.824028 0.0618321\n",
" 0.411746 0.139417\n",
" 0.593904 0.453065\n",
" 0.457076 0.854046\n",
" 0.554123 0.0482456\n",
" 0.434838 0.316422\n",
" 0.899388 0.259067\n",
" 0.0597716 0.516415"
]
},
"execution_count": 117,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"points = rand(Float64, (16,2))\n"
]
},
{
"cell_type": "code",
"execution_count": 118,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 118,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"scatter(points[:,1],points[:,2],xlabel=\"x-coordinate\",ylabel=\"y-coordinate\",legend=false,aspect_ratio =:equal)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Polymake requires a distance matrix to compute the Vietoris rips complex. Here we use the \"Distances\" package to compute a $n \times 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": 119,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"16×16 Array{Float64,2}:\n",
" 0.0 0.670734 1.02139 0.535967 … 0.654219 1.02232 0.345851\n",
" 0.670734 0.0 0.582571 0.330297 0.205125 0.645052 0.358347\n",
" 1.02139 0.582571 0.0 0.487535 0.423885 0.097608 0.845827\n",
" 0.535967 0.330297 0.487535 0.0 0.1666 0.500002 0.393631\n",
" 0.895403 0.986701 0.727044 0.663623 0.789865 0.643458 0.985835\n",
" 0.674019 0.674291 0.526043 0.350196 … 0.481028 0.471475 0.692081\n",
" 0.488818 0.185509 0.69081 0.292566 0.267248 0.733823 0.176352\n",
" 0.560088 0.201051 0.772717 0.418894 0.361351 0.826719 0.214977\n",
" 1.09687 0.588927 0.12257 0.56137 0.465065 0.211141 0.889232\n",
" 0.7977 0.173961 0.427947 0.345074 0.178506 0.502107 0.515765\n",
" 0.661625 0.412742 0.36313 0.144943 … 0.209697 0.361878 0.537876\n",
" 0.383383 0.668668 0.770745 0.371937 0.538083 0.741377 0.521388\n",
" 0.944935 0.342499 0.313982 0.445777 0.29351 0.404541 0.680857\n",
" 0.654219 0.205125 0.423885 0.1666 0.0 0.468078 0.425055\n",
" 1.02232 0.645052 0.097608 0.500002 0.468078 0.0 0.878171\n",
" 0.345851 0.358347 0.845827 0.393631 … 0.425055 0.878171 0.0"
]
},
"execution_count": 119,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"distances = Distances.pairwise(Distances.Euclidean(),points,dims = 1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use Polymake's topaz library to generate the VR complex from our distance matrix. We also have to pass a float, indicating the maximum pairwise distance of vertices, such that they become a face of our abstract simplicial complex. \n",
"Depending on our point set and the value we choose here, we get a higher or lower dimensional abstract simplicial complex. Note that if we'd choose a value of $1$ the generated abstract simplicial complex would contain all faces of the $n-$simplex."
]
},
{
"cell_type": "code",
"execution_count": 130,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"pm::Array>\n",
"{0}\n",
"{1 6 7}\n",
"{1 6 13}\n",
"{1 9 13}\n",
"{2 8 12}\n",
"{2 8 14}\n",
"{3 6 13}\n",
"{3 10 13}\n",
"{4 5}\n",
"{5 10}\n",
"{5 11}\n",
"{6 7 15}\n",
"{9 12 13}\n"
]
},
"execution_count": 130,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vr = topaz.vietoris_rips_complex(distances, 0.33)\n",
"vr.FACETS"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Above we can see the facets, encoded by sets of labeled vertices. Below we see a \n",
"visualization of the VR complex provided by Polymake."
]
},
{
"cell_type": "code",
"execution_count": 150,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Polymake.display_svg(vr)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we use another Polymake function to print the reduced homology of our complex. Each line corresponds to one dimension, starting with dimension $0$. The torsion coefficients are displayed in the curly brackets, which are empty if there is no torsion. The number after the curly brackets is the reduced Betti number."
]
},
{
"cell_type": "code",
"execution_count": 132,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PropertyValue wrapping pm::Array>\n",
"({} 1)\n",
"({} 0)\n",
"({} 0)\n"
]
},
"execution_count": 132,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vr.HOMOLOGY"
]
}
],
"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
}