visualization.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will provide a few examples of multiparameter persistence computations. But first, let's recall the basics! Multiparameter persistence is computed out of **multiparameter persistence modules**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In full generality, a multiparameter persistence module \\\\(M\\\\) is nothing but a family of partially ordered vector spaces, i.e., vector spaces indexed over a poset \\\\(A\\\\), together with morphisms, or linear maps, connecting comparable vector spaces: \n",
    "\n",
    "\\\\[M = \\{\\{M_a\\,:\\,a\\in A\\},\\{\\Phi_{a,b}:M_a\\rightarrow M_b\\}_{a,b\\in A, a\\preceq b}\\}\\\\] \n",
    "\n",
    "such that, for any \\\\(a\\preceq b\\preceq c\\\\), one has \\\\(\\Phi_{a,c}=\\Phi_{b,c}\\circ \\Phi_{a,b}\\\\). \n",
    "\n",
    "Most of the time, these modules are obtained from **multiparameter filtrations**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A multiparameter filtration is the straightforward extension of the usual filtration in the 1D setting. It is a **partially ordered family of spaces**:\n",
    "\n",
    "\\\\[\\mathcal F = \\{X_a\\}_{a\\in A},\\\\]\n",
    "\n",
    "where \\\\( A \\\\) is a poset. The family \\\\(\\mathcal F\\\\) must satisfy \\\\(X_a\\subseteq X_b\\\\) for any \\\\( a,b \\in A\\\\) such that \\\\( a\\preceq b, \\\\), i.e., \\\\(a,b\\\\) are comparable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A very standard example is for \\\\(A=\\mathbb R^d\\\\). Indeed, in that case, there is a natural partial order defined with:\n",
    "\n",
    "\\\\[{\\bf a}\\preceq {\\bf b} \\Longleftrightarrow a_i \\leq b_i,\\ \\ \\ \\forall 1\\leq i \\leq d,\\\\]\n",
    "\n",
    "for any \\\\({\\bf a}, {\\bf b}\\in\\mathbb R^d\\\\). When \\\\(d=2\\\\), one can even represent the multiparameter filtration:\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{@{}c@{\\;}c@{\\;}c@{\\;}c@{\\;}c@{}}\n",
    "X_{1,3} && \\subseteq && X_{2,3} && \\subseteq && X_{3,3}\\\\\n",
    "\\cup | &&  && \\cup | &&  && \\cup |\\\\\n",
    "X_{1,2} && \\subseteq && X_{2,2} && \\subseteq && X_{3,2} \\\\\n",
    "\\cup | &&  && \\cup | &&  && \\cup |\\\\\n",
    "X_{1,1} && \\subseteq && X_{2,1} && \\subseteq && X_{3,1}\n",
    "\\end{array}\n",
    "\\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, it is business as usual, all these inclusions induce morphisms between the homology groups of these spaces:\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{@{}c@{\\;}c@{\\;}c@{\\;}c@{\\;}c@{}}\n",
    "H_*(X_{1,3}) && \\longrightarrow && H_*(X_{2,3}) && \\longrightarrow && H_*(X_{3,3})\\\\\n",
    "\\big\\uparrow &&  && \\big\\uparrow &&  && \\big\\uparrow\\\\\n",
    "H_*(X_{1,2}) && \\longrightarrow && H_*(X_{2,2}) && \\longrightarrow && H_*(X_{3,2}) \\\\\n",
    "\\big\\uparrow &&  && \\big\\uparrow &&  && \\big\\uparrow\\\\\n",
    "H_*(X_{1,1}) && \\longrightarrow && H_*(X_{2,1}) && \\longrightarrow && H_*(X_{3,1})\n",
    "\\end{array}\n",
    "\\\\]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Morever, such multiparameter persistence modules can be compared with the **interleaving distance**:\n",
    "\n",
    "\\\\[d_I(M,M')={\\rm inf}_{\\epsilon >0} M,M'\\text{ are }\\epsilon-\\text{interleaved},\\\\]\n",
    "\n",
    "where \\\\(M,M'\\\\) (with morphisms \\\\(\\Phi,\\Psi\\\\) respectively) are \\\\(\\epsilon\\\\)-interleaved if there are linear maps \\\\(A^{{\\bf a}}_\\epsilon:M_{{\\bf a}}\\rightarrow M_{{\\bf a}+\\epsilon{\\bf 1}}'\\\\) and \\\\(B^{{\\bf a}}_\\epsilon:M'_{{\\bf a}}\\rightarrow M_{{\\bf a}+\\epsilon{\\bf 1}}\\\\) such that \\\\(\\Phi_{{\\bf a},{\\bf a}+2\\epsilon{\\bf 1}}=B^{{\\bf a}+\\epsilon{\\bf 1}}_\\epsilon\\circ A^{{\\bf a}}_\\epsilon\\\\) and \\\\(\\Psi_{{\\bf a},{\\bf a}+2\\epsilon{\\bf 1}}=A^{{\\bf a}+\\epsilon{\\bf 1}}_\\epsilon\\circ B^{{\\bf a}}_\\epsilon\\\\)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next question is obviously: when can we decompose this module into a direct sum of **indecomposable summands** (modules that are not isomorphic to a direct sum of other modules), like in the 1D case? I.e., when do we have:\n",
    "\n",
    "\\\\[H_*(\\mathcal F)=M\\simeq \\oplus_{i=1}^k M_i ??\\\\]\n",
    "\n",
    "Remember, the 1D decomposition was then used to define persistence diagrams, so it would be nice if we could generate an analogue in the multiparameter case. However, the indecomposable summands in the multiparameter setting can be pretty complicated, at least much more than in 1D. For instance, this module is indecomposable:\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{@{}c@{\\;}c@{\\;}c@{\\;}c@{\\;}c@{}}\n",
    "{\\bf k} && [1, 1] && {\\bf k}^2 && [1, 0] && {\\bf k} \\\\\n",
    "[0] &&   && { [0,1]} &&  && [1]\\\\\n",
    "0 && [0] && {\\bf k} && [1] && {\\bf k}\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "and clearly, it has not a simple structure. Actually, it does **not** even contain single copies of the field \\\\( {\\bf k}\\\\) only! Such modules (i.e., modules that contain only single copies \\\\({\\bf k}\\\\)), are called *indicator modules*, and are obviously indecomposable, however they are but just a small subclass of all possible indecomposable modules. Fortunately, we will see in the following that there exists several conditions that can ensure that a given persistence module is decomposable into a direct sum of certain type of indecomposable summands. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, when \\\\(A=\\mathbb Z^d\\\\), the decomposition of a multiparameter persistence module (if it is decomposable) defined on a simplicial complex, can be computed. This is a result of Dey and Xin, presented in [this article](https://arxiv.org/abs/1904.03766#:~:text=Generalized%20Persistence%20Algorithm%20for%20Decomposing%20Multi%2Dparameter%20Persistence%20Modules,-Tamal%20K.&text=Based%20on%20matrix%20reduction%2C%20this,a%201%2Dparameter%20persistence%20module.). The decomposition takes \\\\(O(n^{d(2\\omega+1)}))\\\\) time, where \\\\(n\\\\) is the number of simplices in the simplicial complex, and \\\\(\\omega <2.373\\\\) is the constant for matrix multiplication. This is great, but that complexity can also be pretty large when dealing with simplicial complexes obtained with e.g. Vietoris-Rips filtrations. Also, restricting to \\\\(\\mathbb Z^d\\\\) might also be too restrictive: while it is certainly true that smooth functions, such as Morse functions, only have a finite number of critical values, and thus of homological changes in a multiparameter filtration, it is not necessarily easy to find them beforehand (in order to define the grid in \\\\(\\mathbb Z^d\\\\) that we will use for computing persistence)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we focus on **sublevelset multiparameter persistence modules**, i.e., multiparameter persistence modules coming from two functions defined on the vertices of a simplicial complex \\\\(S\\\\), i.e., we have \\\\(f_1:V(S)\\rightarrow \\mathbb R\\\\) and \\\\(f_2:V(S)\\rightarrow \\mathbb R\\\\). These functions can be (lower-star-)extended to the whole simplicial complex with \\\\(\\tilde f(\\sigma)={\\rm max}\\{f(v)\\,:\\,v\\in\\sigma\\}\\\\). Now, we will study the persistence module induced by the sublevel sets \n",
    "\n",
    "\\\\[M=\\{H_*(\\tilde f_1^{-1}(-\\infty,\\alpha)\\cap \\tilde f_2^{-1}(-\\infty,\\beta))\\}_{\\alpha,\\beta\\in\\mathbb R},\\\\]\n",
    "\n",
    "whose linear maps are induced by the canonical inclusions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our proposed way for getting the decomposition is by slicing the bounding rectangle \\\\( R=[{\\rm min}(\\tilde f_1), {\\rm max}(\\tilde f_1)]\\times [{\\rm min}(\\tilde f_2), {\\rm max}(\\tilde f_2)] \\\\) with a family of lines \\\\(\\mathcal L=\\{\\ell_i\\}_{i=1}^m\\\\). These lines may all originate from the same corner of the rectangle, or have the same slope and be evenly spaced. Each line gives rise to a single-parameter filtration, and has an associated persistence barcode. Then, we can use a black-box matching algorithm to connect the bars of consecutive barcodes together in order to create two-dimensional supports of the summands. Even though this only gives an approximation of the true decomposition of the multiparameter persistence module, the computation time is controlled by the number of lines and can thus be quite fast. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several ways to compute matchings between consecutive barcodes, for instance one can use the matching achieving the bottleneck distance:\n",
    "\n",
    "\\\\[d_B({\\rm dgm}(\\ell_i),{\\rm dgm}(\\ell_{i+1}))= {\\rm inf}_{\\gamma}\\ {\\rm max}_x\\|x-\\gamma(x)\\|_\\infty,\\\\]\n",
    "\n",
    "where \\\\(\\gamma\\\\) ranges over all partial matchings between \\\\({\\rm dgm}(\\ell_i)\\\\) and \\\\({\\rm dgm}(\\ell_{i+1})\\\\). Similarly, one can use the Wasserstein distance:\n",
    "\n",
    "\\\\[d_{W,p}({\\rm dgm}(\\ell_i),{\\rm dgm}(\\ell_{i+1}))^p= {\\rm inf}_{\\gamma}\\ \\sum_x\\|x-\\gamma(x)\\|_\\infty^p.\\\\]\n",
    "\n",
    "A third, more interesting option, is to use vineyards, as described in [this article](https://www.mrzv.org/publications/vineyards/socg06/). Vineyards is just a computation technique that was originally designed to avoid recomputing persistence from scratch when the filtration is slightly changed. More precisely, it allows to update persistence barcodes in linear time (w.r.t. the number of simplices) when only two consecutive simplices are swapped in the filtration. The good news is that it comes with a matching between the associated barcodes (which is basically inferred from the matrix reduction repair after swapping), as a byproduct."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we cythonize the C++ code of the [Dionysus library](https://www.mrzv.org/software/dionysus/examples/pl-vineyard.html#pl-vineyard) that computes vineyards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!cd dionysus_vineyards/ && python setup.py build_ext && find build/ -name \"*.so\" -exec mv {} ../ \\; && rm -r build"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can import the Python package. It depends on the [Gudhi library](http://gudhi.gforge.inria.fr/python/latest/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "from multipers import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are getting close! Let's have a quick look at the code and parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `sublevelsets_multipersistence` computes the decomposition of a given multiparameter filtration. Its most important arguments are:\n",
    "1. `matching`: the black box matching algorithm. You can set it to `'vineyards'` to use the Dionysus library, but any Python function that takes two barcodes as inputs will work fine. For instance, you can use the provided `gudhi_matching` function, which computes the matching based on the 1st Wasserstein distance \\\\(d_{W,1}\\\\).\n",
    "2. `simplextree`: the underlying simplicial complex. It is a SimplexTree from the Gudhi library.\n",
    "3. `filters`: an \\\\(n\\times 2\\\\) NumPy array containing the values of \\\\(f_1\\\\) and \\\\(f_2\\\\). The \\\\(i\\\\)th row of that array should correspond to the filtration values of vertex `[i]` in `simplextree`. \n",
    "4. `homology`: the homology dimension.\n",
    "5. `num_lines`: the number of lines you want to draw in the bounding rectangle \\\\(R\\\\).\n",
    "6. `corner`: the type of family of lines. If `'dg'`: evenly spaced constant slope lines, if `'ur'`: lines originating from upper right corner, if `'ll'`: lines originating from lower left corner.\n",
    "7. `essential`: boolean specifying whether to keep infinite bars (`True`) or not (`False`).\n",
    "8. `bnds_filt`: sides of the bounding rectangle \\\\(R\\\\). If not provided, it is estimated from the data.\n",
    "9. `visu`: boolean used for showing the produced decomposition.\n",
    "10. `plot_per_bar`: boolean used to show each summand individually. Only used if `visu=True`.\n",
    "\n",
    "The other parameters are less important. They specify, e.g., the number of processors when computation are done in parallel, etc. You can run `help(sublevelsets_multipersistence)` for more info. The code then returns, in addition to the decomposition, the family of lines and the bounding rectangle (that was potentially estimated from the data)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two modules with the same rank invariant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A very standard tool to analyze multiparameter persistence modules is the so-called **rank invariant**. It is defined as \n",
    "\n",
    "\\\\[R(M):\\left\\{\\begin{align*}\n",
    "\\Omega & \\rightarrow\\mathbb N \\\\\n",
    "(a,b) & \\mapsto {\\rm rank}(M_a \\rightarrow M_b)={\\rm dim}({\\rm im}(M_a \\rightarrow M_b))\n",
    "\\end{align*}\n",
    "\\right.\\\\]\n",
    "where \\\\(\\Omega=\\{(a,b)\\in A\\times A\\,:\\, a\\preceq b\\}\\\\)\n",
    "\n",
    "However, it is known that it does not contain more information than the family of barcodes (i.e., it ignores the matching between consecutive barcodes), so it is quite easy to generate different persistence modules with the same rank invariant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an extremely simple simplicial complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  0 --- 3 / 2\n",
    "#        |\n",
    "#        |\n",
    "#        1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplextree = gd.SimplexTree()\n",
    "simplextree.insert([0,3])\n",
    "simplextree.insert([1,3])\n",
    "simplextree.insert([2])\n",
    "F1 = np.array([0,1,1,1])\n",
    "F2 = np.array([1,0,1,1])\n",
    "filters = np.hstack([F1[:,np.newaxis], F2[:,np.newaxis]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see its decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "D, A, B, _ = sublevelsets_multipersistence(\n",
    "    \"vineyards\", simplextree, filters, homology=0, num_lines=100, essential=True, corner='dg',\n",
    "    bnds_filt=[0,2,0,2], visu=True, plot_per_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, on the diagonal, the two bars in the associated persistence barcode are exactly the same. Hence, if we do not use a matching function between barcodes that uses prior knowledge about how the barcodes were generated (as vineyards do), the matching will necessarily be arbitrary, and potentially lead to wrong decompositions. Here is, for instance, the matching produced by the 1st Wasserstein distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D, A, B, _ = sublevelsets_multipersistence(\n",
    "    gudhi_matching, simplextree, filters, homology=0, num_lines=100, essential=True, corner='dg',\n",
    "    bnds_filt=[0,2,0,2], visu=True, plot_per_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks definitely wrong. Vineyards can really make a difference!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can go one step further and compute a vectorization of this decomposition. We want something more powerful than the rank invariant so we need to take the colors of the bars into account somehow. We can do this simply by centering Gaussian function on all lines, and weight them by the area of the subsets produced by the family of lines. \n",
    "\n",
    "\\\\[I_{\\mathcal D,\\{\\ell_i\\}_i}(P)=\\sum_{S\\in \\mathcal D}w_s(S)\\left(w_l(\\ell^*)\\cdot{\\rm min}_{i}{\\rm exp}\\left(-\\frac{d(P,\\ell_i)^2}{\\sigma^2}\\right)\\right)\\\\]\n",
    "\n",
    "This is something you can do with the function `multipersistence_image`. The important parameters are:\n",
    "\n",
    "1. `decomposition`: the decomposition as produced by `sublevelsets_multipersistence` or `interlevelsets_multipersistence`.\n",
    "2. `bnds`: the sides of the bounding rectange \\\\(R\\\\).\n",
    "3. `resolution`: the image resolution, i.e., the number of pixels on each side of the image.\n",
    "4. `bandwidth`: the bandwidth of the Gaussian functions.\n",
    "5. `power`: the Gaussian functions are weighted by the summand area, raised to that power: \\\\(w_s(S)={\\rm area}(S)^p\\\\)\n",
    "6. `line_weight`: a function that weights the Gaussian functions depending on the slopes of the lines. Any function that takes two numbers as input (the coordinates of a unit directing vector) works fine. This is useful only if the decomposition was computed with `corner='ur'` or `corner='ll'`.\n",
    "\n",
    "Weighting by the summand area ensures more discriminativity than the rank invariant, and weighting by the line slopes ensures stability. Indeed, a results of Landi in [this article](https://arxiv.org/pdf/1412.3374.pdf) shows that the barcodes associated to lines drawn on a multiparameter persistence module are upper bounded by the interleaving distance times the inverse of the line slope (which makes sense since the interleaving distance only measures differences along lines with slope 1):\n",
    "\n",
    "\\\\[d_B({\\rm dgm}^M(\\ell),{\\rm dgm}^{M'}(\\ell)) \\leq \\frac{1}{{\\min(e^\\ell_0,e^\\ell_1)}} d_I(M,M'),\\\\]\n",
    "\n",
    "where \\\\({\\bf e^\\ell}=[e^\\ell_0,e^\\ell_1]\\\\) is a unit directing vector for \\\\(\\ell\\\\)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = multipersistence_image(D, np.array(B), resolution=[200,200], bandwidth=.1, \n",
    "                           power=1., line_weight=lambda x: min(x[0],x[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.flip(I, 0), vmin=0., vmax=2.)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also provide code to compute the [extension of landscapes](https://jmlr.org/papers/v21/19-054.html) to multiparameter persistence modules. They are simply computed by aggregating all the 1D persistence landscapes along all lines. The function `multiparameter_landscape` can compute such an extension when the decomposition was computed with `corner='dg'`. Its important parameters are:\n",
    "\n",
    "1. `decomposition`: the decomposition as produced by `sublevelsets_multipersistence` or `interlevelsets_multipersistence`.\n",
    "2. `bnds`: the sides of bounding rectange \\\\(R\\\\).\n",
    "3. `delta`: the distance between consecutive lines. It can be computed from the second output of `sublevelset_multipersistence`.\n",
    "4. `resolution`: number of pixels for each landscape axis.\n",
    "5. `k`: the number of landscapes to compute. If `None`, returns the sum of all tent functions, i.e., the [silhouette](https://geometrica.saclay.inria.fr/team/Fred.Chazal/papers/cflrw-scpls-14/cflrw-scpls-14.pdf).\n",
    "6. `power`: the tent functions are raised to that power. This is useful only if `k=None`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = np.abs(A[0,0] - A[1,0]) if A[0,0] != A[1,0] else np.abs(A[2,2] - A[1,2])\n",
    "L = multipersistence_landscape(D, B, delta, resolution=[200,200], k=None, power=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "#plt.imshow(np.flip(L.sum(axis=2), 0), vmin=0., vmax=3.)\n",
    "plt.imshow(np.flip(L, 0), vmin=0., vmax=3.)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's check that rank invariants and multiparameter persistence landscapes are indeed less discriminative than multiparameter persistence images with this second simple simplicial complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  0\n",
    "#\n",
    "#\n",
    "#        1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplextree = gd.SimplexTree()\n",
    "simplextree.insert([0])\n",
    "simplextree.insert([1])\n",
    "F1 = np.array([0,1])[:,np.newaxis]\n",
    "F2 = np.array([1,0])[:,np.newaxis]\n",
    "filters = np.hstack([F1,F2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The decomposition is not the same since the multiparameter persistence modules are different. However the barcodes associated to the lines are all the same, which means that these modules are the same from the point of view of the rank invariant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "D, A, B, _ = sublevelsets_multipersistence(\n",
    "    'vineyards', simplextree, filters, homology=0, num_lines=100, essential=True, \n",
    "    bnds_filt=[0,2,0,2], epsilon=1e-10, visu=True, plot_per_bar=False, bnds_visu=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute our two descriptors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = np.abs(A[0,0] - A[1,0]) if A[0,0] != A[1,0] else np.abs(A[2,2] - A[1,2])\n",
    "I = multipersistence_image(D, B, resolution=[200,200], bandwidth=.1, power=0.)\n",
    "L = multipersistence_landscape(D, B, delta, resolution=[200,200], k=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is clear that the multiparameter persistence images are different while the multiparameter persistence landscapes are the same!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.flip(I, 0), vmin=0., vmax=2.)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.flip(L.sum(axis=2), 0), vmin=0., vmax=3.)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Indecomposable persistence module"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is great, but what happens when the code is applied on a multiparameter persistence module that is not indecomposable into indicators? Then, the code will produce a decomposition into a direct sum of indicators, which must then be wrong. Indeed, let's look at the following simplicial complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#         5\n",
    "#       / | \\\n",
    "#      /  4  \\\n",
    "#     / / | \\ \\\n",
    "#   1 --- 2 --- 3\n",
    "#     \\   |   /\n",
    "#      \\  |  /\n",
    "#       \\ | /\n",
    "#         0     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplextree = gd.SimplexTree()\n",
    "simplextree.insert([0])\n",
    "simplextree.insert([1])\n",
    "simplextree.insert([2])\n",
    "simplextree.insert([3])\n",
    "simplextree.insert([4])\n",
    "simplextree.insert([5])\n",
    "simplextree.insert([0, 1])\n",
    "simplextree.insert([0, 3])\n",
    "simplextree.insert([1, 2])\n",
    "simplextree.insert([2, 3])\n",
    "simplextree.insert([2, 4])\n",
    "simplextree.insert([1, 4])\n",
    "simplextree.insert([3, 4])\n",
    "simplextree.insert([1, 5])\n",
    "simplextree.insert([3, 5])\n",
    "simplextree.insert([4, 5])\n",
    "simplextree.insert([1, 2, 4])\n",
    "simplextree.insert([2, 3, 4])\n",
    "simplextree.insert([1, 4, 5])\n",
    "simplextree.insert([3, 4, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This complex filtered by the following filtration values has the multiparameter persistence module (in homology dimension 1) that was described at the beginning of the notebook. \n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{@{}c@{\\;}c@{\\;}c@{\\;}c@{\\;}c@{}}\n",
    "{\\bf k} && [1, 1] && {\\bf k}^2 && [1, 0] && {\\bf k} \\\\\n",
    "[0] &&   && { [0,1]} &&  && [1]\\\\\n",
    "0 && [0] && {\\bf k} && [1] && {\\bf k}\n",
    "\\end{array}\n",
    "\\\\]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtration = np.array([[1, 1],[1, 2],[2, 2],[1, 2],[3, 2.5],[1, 3]])\n",
    "print(filtration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still, we can compute a decomposition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "D, A, B, _ = sublevelsets_multipersistence(\n",
    "    \"vineyards\", simplextree, filtration, homology=1, num_lines=100, essential=True, \n",
    "    corner=\"dg\", bnds_filt=[0,4,1.5,4], visu=True, plot_per_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute associated descriptors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = np.abs(A[0,0] - A[1,0]) if A[0,0] != A[1,0] else np.abs(A[2,2] - A[1,2])\n",
    "I = multipersistence_image(D, np.array(B), resolution=[200,200], bandwidth=.1)\n",
    "L = multipersistence_landscape(D, B, delta, resolution=[200,200], k=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.flip(I, 0))\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.flip(L[:,:,0], 0))\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.imshow(np.flip(L[:,:,1], 0))\n",
    "plt.show()\n",
    "plt.figure()\n",
    "plt.imshow(np.flip(L[:,:,2], 0))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The message here is that one has to be careful when using this code if one has no prior about the decomposability of its data, otherwise the code will produce decomposition that do not necessarily correspond to the data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Interlevelsets multipersistence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A next natural question is then: when can we ensure that multiparameter persistence modules are decomposable? It turns out to be a tough question. For instance, in [this article](https://arxiv.org/abs/1812.05261), the authors provide a way to check whether a given multiparameter persistence module is decomposable into a specific class of indicator summands called *interval modules*, and in [this article](https://arxiv.org/abs/2002.08894), the authors characterize the class of modules that are decomposable into *rectangle indicator modules*.\n",
    "\n",
    "Another very important setting is the case of **exact** multiparameter persistence modules, as presented in [this article](https://link.springer.com/article/10.1007/s00454-019-00165-z). Being exact simply means that the sequence \\\\(M_{x,y}\\stackrel{\\phi}{\\rightarrow} M_{x',y}\\oplus M_{x,y'}\\stackrel{\\psi}{\\rightarrow} M_{x',y'}\\\\) is exact, i.e., \\\\({\\rm im}(\\phi)\\simeq{\\rm ker}(\\psi)\\\\), where \\\\(\\phi=\\Phi_{(x,y),(x',y)}\\oplus\\Phi_{(x,y),(x,y')}\\\\) and \\\\(\\psi=\\Phi_{(x',y'),(x',y)}-\\Phi_{(x',y'),(x,y')}\\\\). Practically speaking, this ensures that the following diagram commutes:\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{@{}c@{\\;}c@{\\;}c@{\\;}c@{\\;}c@{}}\n",
    "M_{x,y'} && \\longrightarrow && M_{x',y'}  \\\\\n",
    "\\big\\uparrow &&  && \\big\\uparrow \\\\\n",
    "M_{x,y} && \\longrightarrow && M_{x',y}\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "and that each element in \\\\(M_{(x',y')}\\\\) that has preimages in both \\\\(M_{(x',y)}\\\\) and \\\\(M_{(x,y')}\\\\), has a preimage common to both in \\\\(M_{(x,y)}\\\\)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A very nice result about exact persistence bimodules is that they can always be decomposed into **block indicator summands**, i.e., indicator summands of one of the four following forms:\n",
    "\n",
    "<img src=\"images/blocks.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But when exactly is a multiparameter persistence module exact? It can actually be shown that this happens for **interlevelset multiparameter persistence modules**:\n",
    "\n",
    "\\\\[M=\\{H_*(\\tilde f^{-1}([\\alpha,\\beta])\\}_{\\alpha,\\beta\\in\\mathbb R, \\alpha\\leq\\beta},\\\\]\n",
    "\n",
    "Let's check this on an example! Here is a simplicial complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#      5\n",
    "#      |\n",
    "#      4   7\n",
    "#     / \\ /\n",
    "#    2   3\n",
    "#   / \\ /\n",
    "#  6   1\n",
    "#      |\n",
    "#      0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplextree = gd.SimplexTree()\n",
    "simplextree.insert([0])\n",
    "simplextree.insert([1])\n",
    "simplextree.insert([2])\n",
    "simplextree.insert([3])\n",
    "simplextree.insert([4])\n",
    "simplextree.insert([5])\n",
    "simplextree.insert([6])\n",
    "simplextree.insert([7])\n",
    "simplextree.insert([0, 1])\n",
    "simplextree.insert([1, 2])\n",
    "simplextree.insert([1, 3])\n",
    "simplextree.insert([2, 4])\n",
    "simplextree.insert([3, 4])\n",
    "simplextree.insert([4, 5])\n",
    "simplextree.insert([2, 6])\n",
    "simplextree.insert([3, 7])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's define a vertex-valued function (we only need one now for interlevelsets)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtration = np.array([0,1,2,2,3,4,0.4,3.6])\n",
    "print(filtration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Decomposition of interlevelset multiparameter persistence modules can be done with the `interlevelsets_multipersistence` function. Its parameters are roughly the same than `sublevelsets_multipersistence` except that, since it is only defined on \\\\(\\{(x,y)\\,:\\,y\\geq x\\}\\\\) it also requires a `basepoint`, i.e., a point on the diagonal from which all lines will originate. Also, for practical and esthetic reasons, we actually provide the symmetric (across the \\\\(y\\\\)-axis) of the decomposition. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D0, A, B, _ = interlevelsets_multipersistence(\n",
    "    \"vineyards\", simplextree, filtration, homology=0, num_lines=100, basepoint=2.5, essential=True,\n",
    "    visu=True, plot_per_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mmmmh... those summands are similar, but not quite equal to blocks. What's going on?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason why we do not quite have blocks is because we made an approximation by extending the vertex-valued function to all simplices, remember? Hence, the module is not exactly an interlevelset multiparameter persistence module. To check that this is indeed what's happening, we provide a function `barycentric_subdivision`, that can refine a simplicial complex in the following way:\n",
    "\n",
    "<img src=\"images/subdiv.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsub = 4\n",
    "singlefiltration = filtration\n",
    "subdivision = gd.SimplexTree()\n",
    "[subdivision.insert(s,f) for s,f in simplextree.get_filtration()]\n",
    "for _ in range(nsub): \n",
    "    for s,_ in subdivision.get_filtration():\n",
    "        if len(s) == 1:\n",
    "            subdivision.assign_filtration(s,singlefiltration[s[0]])\n",
    "        else:\n",
    "            subdivision.assign_filtration(s,-1e10)\n",
    "    subdivision.make_filtration_non_decreasing()\n",
    "    list_splx = list(subdivision.get_filtration())\n",
    "\n",
    "    spt = barycentric_subdivision(subdivision, list_splx, use_mean=True)\n",
    "    ftr = np.array([spt.filtration(s) for s,_ in spt.get_skeleton(0)])\n",
    "\n",
    "    subdivision = spt\n",
    "    singlefiltration = ftr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D0, A, B, _ = interlevelsets_multipersistence(\n",
    "    \"vineyards\", subdivision, singlefiltration, homology=0, num_lines=100, basepoint=2.5, essential=True,\n",
    "     visu=True, plot_per_bar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yep, now we have blocks! The math actually works!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our last example is very common in TDA, and deals with the limitation of geometric complexes such as Vietoris-Rips or Alpha. Indeed, even though it was shown in [this article](https://arxiv.org/abs/1207.3885) that they are stable w.r.t. the Gromov-Hausdorff distance, this distance is itself very sensitive to outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's process the following point cloud, generated from a noisy circle plus three outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npts = 100\n",
    "data = np.random.uniform(low=0, high=2*np.pi, size=[npts,1])\n",
    "noise = np.random.uniform(low=-.3, high=.3, size=[npts,2])\n",
    "X = np.hstack([np.cos(data),np.sin(data)]) + noise\n",
    "X = np.vstack([X, np.array([[-.1,.2],[.3,.1],[-.2,-.1]])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.scatter(X[:,0],X[:,1])\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because of these outliers, the Alpha persistence diagram is actually **not** stable, in the sense that these outliers fool persistent homology into thinking that there are more than just one cycle. Indeed, there are more than one point that are far away from the diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_complex = gd.AlphaComplex(points=X)\n",
    "simplex_tree = alpha_complex.create_simplex_tree(max_alpha_square=2)\n",
    "simplex_tree.persistence()\n",
    "dgm_alpha = simplex_tree.persistence_intervals_in_dimension(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.scatter(dgm_alpha[:,0], dgm_alpha[:,1])\n",
    "plt.plot([dgm_alpha.min(), dgm_alpha.max()],[dgm_alpha.min(), dgm_alpha.max()])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Howver, these outliers can be detected with simple tools like kernel density estimation since the density estimates are quite small w.r.t. the other points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kde = KernelDensity(kernel='gaussian', bandwidth=.3).fit(X)\n",
    "density = kde.score_samples(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.scatter(X[:,0],X[:,1], c=density)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hence, why not filtering jointly by geometry **and** density? There is one limitation though: Alpha filtrations are not sublevelset filtrations of some lower-star function. Fortunately, every filtration can be turned into the sublevelset filtration of a function using, again, barycentric subdivision. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_splx = simplex_tree.get_filtration()\n",
    "subdivision = barycentric_subdivision(simplex_tree)\n",
    "filtration_rips = np.array([subdivision.filtration(s) for s,_ in subdivision.get_skeleton(0)])[:,None]\n",
    "for s,_ in simplex_tree.get_filtration():\n",
    "    if len(s) == 1:\n",
    "        simplex_tree.assign_filtration(s, -density[s[0]])\n",
    "    else:\n",
    "        simplex_tree.assign_filtration(s, -1e10)\n",
    "simplex_tree.make_filtration_non_decreasing()\n",
    "subdivision = barycentric_subdivision(simplex_tree, list_splx)\n",
    "filtration_dens = np.array([subdivision.filtration(s) for s,_ in subdivision.get_skeleton(0)])[:,None]\n",
    "filtration = np.hstack([filtration_rips, filtration_dens])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready. Let's check that decomposition!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D1, A, B, _ = sublevelsets_multipersistence(\n",
    "    \"vineyards\", subdivision, filtration, homology=1, num_lines=300, essential=False,\n",
    "    bnds_filt=[min(filtration[:,0]),max(filtration[:,0]),min(filtration[:,1]),max(filtration[:,1])], \n",
    "    visu=True, plot_per_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is definitely one big summand for the circle and only small summands for the rest of the decomposition. This indicates that this kind of construction is robust, which was recently proved true in [this article](https://arxiv.org/abs/2010.09628)."
   ]
  },
  {
   "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.9.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}