{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MVA 2021-22  - TP 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To download this notebook:\n",
    "\n",
    "http://geometrica.saclay.inria.fr/team/Fred.Chazal/MVA2022.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can install gudhi using pip or conda:\n",
    "\n",
    "https://gudhi.inria.fr/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gudhi as gd\n",
    "print(gd.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this first TP is to get you familiar with the basic data structures in GUDHI to  build and manipulate simplicial complexes and filtrations. "
   ]
  },
  {
   "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\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.cm\n",
    "\n",
    "#%matplotlib inline\n",
    "%matplotlib notebook"
   ]
  },
  {
   "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",
    "L = st.get_filtration()\n",
    "for simplex in L:\n",
    "    print(simplex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filtration value of a simplex can be changed in the following way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st.assign_filtration((2,3,4),1.0)\n",
    "\n",
    "L = st.get_filtration()\n",
    "for simplex in L:\n",
    "    print(simplex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Warning! ** Take care that after changing the filtration value of a simplex, the result could no longer be a filtration, as illustrated below :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Giving the edge [3,4] the value 1.5:\")\n",
    "st.assign_filtration((3,4),1.5)\n",
    "L = st.get_filtration()\n",
    "for simplex in L:\n",
    "    print(simplex)\n",
    "print(\"The result is no longer a filtration : [3,4] has a higher value than its coface [2,3,4]\")\n",
    "print(\"To fix the problem, use make_filtration_non_decreasing()\")\n",
    "st.make_filtration_non_decreasing()\n",
    "L = st.get_filtration()\n",
    "for simplex in L:\n",
    "    print(simplex)"
   ]
  },
  {
   "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 value of [1,2]=\",st.filtration([1,2]))\n",
    "print(\"filtration value of [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))\n",
    "L = st.get_skeleton(1)\n",
    "for simplex in L:\n",
    "    print(simplex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1. \n",
    "Make a few experiments with the simplex tree functions ( https://gudhi.inria.fr/python/latest/simplex_tree_ref.html ), e.g. changing the filtrations values, trying to assign values to simplices that do not lead to a filtration,... And observe the effects on the filtration. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtrations, persistence and Betti numbers computation"
   ]
  },
  {
   "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",
    "L = st.get_filtration()\n",
    "for simplex in L:\n",
    "    print(simplex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before computing Betti numbers, we first need to compute persistence of the 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",
    "## Here, for the moment, we use it as a preprocessing step to compute Betti numbers. \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": [
    "print(st.betti_numbers())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2. \n",
    "Define another filtration of the simplicial complex and check that the choice of the filtration does not change the betti numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3. Persistence for functions defined over a grid\n",
    "When a function is defined over a grid on $[0,1]$, $[0,1]^2$, ... the grid can be directly used to build the filtration using a cubical complex as illustrated in the following examples. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Persistence of 1D function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let consider $f: t \\mapsto sin(2t)+sin(3t)$ defined over $[0, 2\\pi]$. \n",
    "\n",
    "Build a table with 200 values of f between 0 and $2\\pi$. Plot the function, compute the persistence diagram of its sublevelsets, and draw its persistence diagram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Remark:\n",
    "The name top_dimensional_cells is because gudhi gives the grid values to top-dimensional cells and deduces values for other cells, instead of giving values to vertices and deducing values for other cells. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Persistence of 2D function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $p_0=(0.25, 0.25), p_1=(0.75, 0.75), p_2 = (0.0, 1.0)$ and $p_3 = (1.0, 0.0)$ be 4 points in the plane $\\mathbb{R}^2$ and $\\sigma=0.05$.\n",
    "1. Build on such a complex the sublevelset filtration of the function \n",
    "$$f(p)=\\exp(-\\frac{\\|p-p_0\\|^2}{\\sigma})+3\\exp(-\\frac{\\|p-p_1\\|^2}{\\sigma}) - 4*\\exp(-\\frac{\\|p-p_2\\|^2}{\\sigma}) \n",
    "- 2 \\exp(-\\frac{\\|p-p_3\\|^2}{\\sigma})$$ \n",
    "defined over $[-0.5,1.5]^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the persistence diagram of the sublevel set filtration of $f$ and compute the persistence diagram of the upperlevel set filtration of $f$ and compare the obtained diagram to the previous one. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Persistence over a filtered simplicial complex and Betti number. \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.\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": "markdown",
   "metadata": {},
   "source": [
    "# Vietoris-Rips and alpha-complex filtrations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the definition of Vietoris-Rips and $\\alpha$-complexes, see the slides of the course: https://geometrica.saclay.inria.fr/team/Fred.Chazal/slides/Persistence2022.pdf\n",
    "\n",
    "See also the following book, p.137\n",
    "https://hal.inria.fr/hal-01615863v2/document\n",
    "\n",
    "Take care that in GUDHI the α-complex filtration is indexed by the square of the radius of the smallest empty circumscribing ball. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are basic instructions to build Vietoris-Rips and α-complex filtrations (and compute their persistent homology)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Random point cloud"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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",
    "dgm = simplex_tree.persistence(homology_coeff_field=2, min_persistence=0)\n",
    "gd.plot_persistence_diagram(dgm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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 4.  Sampled torus in $\\mathbb{R}^4$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Randomly sample $n$ points (try different values of $n$) on the paramatrized torus in $\\mathbb{R}^4$:\n",
    "\n",
    "$$(s,t) \\to (\\cos(s),\\sin(s),\\cos(t),\\sin(t)), \\ \\  \\ (s,t) \\in [0,2\\pi] \\times [0,2\\pi]$$\n",
    "\n",
    "and compute the persistence diagrams of the resulting $\\alpha$-complex filtration. \n",
    "\n",
    "+ Do the same with the Vietoris-Rips complex. \n",
    "\n",
    "+ What do you observe?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Now, sample the points along a 1D line embedded in the torus according to the following parametrization:\n",
    "\n",
    "$$t \\to (\\cos(t),\\sin(t),\\cos(5t),\\sin(5t)), \\ \\  \\ t \\in [0,2\\pi] $$\n",
    "\n",
    "and do the same experiment as previously. What do you observe? Explain it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5. \n",
    "+ Randomly sample n = 100 points on the unit circle in the Euclidean plane.\n",
    "+ For R in np.arange(0.0,0.5,0.01), compute the Betti numbers of the subcomplex of the Rips-Vietoris filtration (up to dimension 2) made of the simplices with index value at most R and plot the curve giving the Betti numbers as functions of R. These curves are called the Betti curves of the filtration. \n",
    "+ Can we get the same curves directly from the persistence diagram of the Rips-Vietoris filtration (you will have to guess what the persistence diagrams represent)? If so, compute them using the persistence diagram. \n",
    "+ Same questions using the α-complex filtrations (find a right range of values for α), and try to increase the number of points in the initial point cloud. \n",
    "+ Do the same for the point cloud sampled on the 2D torus in $\\mathbb{R}^4$ from the above exercise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.13"
  },
  "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
}