{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MVA 2018-19"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To download this notebook or its pdf version:\n",
    "\n",
    "http://geometrica.saclay.inria.fr/team/Fred.Chazal/MVA2018.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Documentation for the latest version of Gudhi: \n",
    "\n",
    "http://gudhi.gforge.inria.fr/python/latest/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simplicial complexes and simplex trees"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Gudhi, (filtered) simplicial complexes are encoded through a data structure called simplex tree. Here is a very simple example illustrating the use of simplex tree to represent simplicial complexes. See the Gudhi documentation for a complete list of functionalities. Try the following code and a few other functionalities from the documentation to get used to the Simplex Tree data structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import gudhi as gd\n",
    "import random as rd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st = gd.SimplexTree() # Create an empty simplicial complex\n",
    "\n",
    "# Simplicies can be inserted 1 by 1\n",
    "# Vertices are indexed by integers\n",
    "if st.insert([0,1]):\n",
    "    print(\"First simplex inserted!\")\n",
    "st.insert([1,2])\n",
    "st.insert([2,3])\n",
    "st.insert([3,0])\n",
    "st.insert([0,2])\n",
    "st.insert([3,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = st.get_filtration() # Get a list with all simplices\n",
    "# Notice that inserting an edge automatically inserts its vertices, if they were not already in the complex\n",
    "for simplex in L:\n",
    "    print(simplex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Insert the 2-skeleton, giving some filtration values to the faces\n",
    "st.insert([0,1,2],filtration=0.1)\n",
    "st.insert([1,2,3],filtration=0.2)\n",
    "st.insert([0,2,3],filtration=0.3)\n",
    "st.insert([0,1,3],filtration=0.4)\n",
    "\n",
    "# If you add a new simplex with a given filtration value, all its faces that \n",
    "# were not in the complex are inserted with the same filtration value\n",
    "st.insert([2,3,4],filtration=0.7)\n",
    "st.get_filtration()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Many operations can be done on simplicial complexes, see also the Gudhi documentation and examples\n",
    "print(\"dimension=\",st.dimension())\n",
    "print(\"filtration[1,2]=\",st.filtration([1,2]))\n",
    "print(\"filtration[4,2]=\",st.filtration([4,2]))\n",
    "print(\"num_simplices=\", st.num_simplices())\n",
    "print(\"num_vertices=\", st.num_vertices())\n",
    "print(\"skeleton[2]=\", st.get_skeleton(2))\n",
    "print(\"skeleton[1]=\", st.get_skeleton(1))\n",
    "print(\"skeleton[0]=\", st.get_skeleton(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtrations and persistence diagrams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As an example, we assign to each simplex its dimension as filtration value\n",
    "for splx in st.get_filtration():\n",
    "    st.assign_filtration(splx[0],len(splx[0])-1)\n",
    "# Let the structure know that we have messed with it and an old filtration cache may be invalid.\n",
    "# This is redundant here because get_filtration() does it anyway, but not all functions do.\n",
    "st.initialize_filtration()\n",
    "st.get_filtration()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To compute the persistence diagram of the filtered complex\n",
    "# By default it stops at dimension-1, use persistence_dim_max=True\n",
    "# to compute homology in all dimensions\n",
    "diag = st.persistence(persistence_dim_max=True)\n",
    "# Display each interval as (dimension, (birth, death))\n",
    "print(diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot a persistence diagram\n",
    "gd.plot_persistence_diagram(diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gd.plot_persistence_barcode(diag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To compute the bottleneck distance between diagrams\n",
    "st2 = gd.SimplexTree()\n",
    "st2.insert([0,1],filtration=0.0)\n",
    "st2.insert([1,2],filtration=0.1)\n",
    "st2.insert([2,0],filtration=0.2)\n",
    "st2.insert([0,1,2],filtration=0.5)\n",
    "diag2 = st2.persistence()\n",
    "gd.plot_persistence_diagram(diag2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diag_0 = st.persistence_intervals_in_dimension(0)\n",
    "diag2_0 = st2.persistence_intervals_in_dimension(0)\n",
    "dB0 = gd.bottleneck_distance(diag_0, diag2_0)\n",
    "\n",
    "diag_1 = st.persistence_intervals_in_dimension(1)\n",
    "diag2_1 = st2.persistence_intervals_in_dimension(1)\n",
    "dB1 = gd.bottleneck_distance(diag_1, diag2_1)\n",
    "\n",
    "print(\"bottleneck distance in dimension 0:\", dB0)\n",
    "print(\"bottleneck distance in dimension 1:\", dB1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "1. Recall the torus is homeomorphic to the surface obtained by identifying the opposite sides of a square as illustrated below. ![Figure 1](TorusTriangle.png) Using Gudhi, construct a triangulation (2-dimensional simplicial complex) of the Torus. Define a filtration on it, compute its persistence and use it to deduce the Betti numbers of the torus (check that you get the correct result using the function betti_numbers()).\n",
    "2. Use Gudhi to compute the Betti numbers of a sphere of dimension 2 and of a sphere of dimension 3 (hint: the k -dimensional sphere is homeomorphic to the boundary of a (k+1)-dimensional simplex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stability of persistence for functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this exercise is to illustrate the persistence stability theorem for functions on a very simple example. The code below allows to define a simplicial complex (the so-called α-complex) triangulating a set of random points in the unit square in the plane. Although we are not using it for this course, note that gudhi also provides the grid-like CubicalComplex, which may be a more natural choice to represent a function on a square."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pts=1000\n",
    "# Build a random set of points in the unit square\n",
    "X = np.random.rand(n_pts, 2)\n",
    "# Compute the alpha-complex filtration\n",
    "alpha_complex = gd.AlphaComplex(points=X)\n",
    "st = alpha_complex.create_simplex_tree()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $p_0=(0.25, 0.25)$ and $p_1=(0.75, 0.75)$ be 2 points in the plane $\\mathbb{R}^2$ and $\\sigma=0.05$.\n",
    "1. Build on such a complex the sublevelset filtration of the function $$f(p)=\\exp(-\\frac{\\|p-p_0\\|^2}{\\sigma})+3\\exp(-\\frac{\\|p-p_1\\|^2}{\\sigma})$$ and compute its persistence diagrams in dimensions 0 and 1.\n",
    "2. Compute the persistence diagrams of random perturbations of f and compute the Bottleneck distance between these persistence diagrams and the perturbated ones. Verify that the persistence stability theorem for functions is satisfied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vietoris-Rips and alpha-complex filtrations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are basic instructions to build Vietoris-Rips and α-complex filtrations and compute their persistence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a random point cloud in 3D\n",
    "nb_pts=100\n",
    "pt_cloud = np.random.rand(nb_pts,3)\n",
    "#Build Rips-Vietoris filtration and compute its persistence diagram\n",
    "rips_complex = gd.RipsComplex(pt_cloud,max_edge_length=0.5)\n",
    "simplex_tree = rips_complex.create_simplex_tree(max_dimension=3)\n",
    "print(\"Number of simplices in the V-R complex: \",simplex_tree.num_simplices())\n",
    "diag = simplex_tree.persistence(homology_coeff_field=2, min_persistence=0)\n",
    "gd.plot_persistence_diagram(diag).show()\n",
    "#Compute Rips-Vietoris filtration and compute its persistence diagram from\n",
    "#a pairwise distance matrix\n",
    "dist_mat = []\n",
    "for i in range(nb_pts):\n",
    "    ld = []\n",
    "    for j in range(i):\n",
    "        ld.append(np.linalg.norm(pt_cloud[i,:]-pt_cloud[j,:]))\n",
    "    dist_mat.append(ld)\n",
    "rips_complex2 = gd.RipsComplex(distance_matrix=dist_mat,max_edge_length=0.5)\n",
    "simplex_tree2 = rips_complex2.create_simplex_tree(max_dimension=3)\n",
    "diag2 = simplex_tree2.persistence(homology_coeff_field=2, min_persistence=0)\n",
    "gd.plot_persistence_diagram(diag2).show()\n",
    "#Compute the alpha-complex filtration and compute its persistence\n",
    "alpha_complex = gd.AlphaComplex(points=pt_cloud)\n",
    "simplex_tree3 = alpha_complex.create_simplex_tree()\n",
    "print(\"Number of simplices in the alpha-complex: \",simplex_tree3.num_simplices())\n",
    "diag3 = simplex_tree3.persistence(homology_coeff_field=2, min_persistence=0)\n",
    "gd.plot_persistence_diagram(diag3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Illustrate the stability theorem for persistence diagrams of Vietoris-Rips and α-complex filtrations (take into accound that AlphaComplex uses the square of distances as filtration values).\n",
    "2. What happens to Vietoris-Rips and α-complex filtrations when the size of the point cloud increases? When the ambient dimension increases?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}