{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5e027782",
   "metadata": {},
   "source": [
    "# TDA: TP 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba984f4b",
   "metadata": {},
   "source": [
    "For these exercises, we will use the library [Gudhi](https://gudhi.inria.fr/). You can install it with `conda` or `pip` (on new Apple laptops, only `conda` is available). Let us check that you have a recent version (current is 3.6.0). If for some reason you end up with an old version, you can try to install `gudhi=3.6.0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "976e38fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gudhi as gd\n",
    "print(gd.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89ea55eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will need a few other libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.cm\n",
    "# To get nice interactive plots. May require installing ipympl and restarting the kernel.\n",
    "# Other possibilities than notebook include qt, osx, tk... If nothing works, just remove this line\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e783b8c",
   "metadata": {},
   "source": [
    "The documentation for the Python interface of Gudhi is [here](https://gudhi.inria.fr/python/latest/). There are also some [tutorials](https://github.com/GUDHI/TDA-tutorial), an [issue tracker](https://github.com/GUDHI/gudhi-devel/issues), a [mailing-list](https://sympa.inria.fr/sympa/arc/gudhi-users/), etc."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49e4cdd1",
   "metadata": {},
   "source": [
    "## Functions, cubical complex\n",
    "### Volcano\n",
    "Let us first define a function from $\\mathbb{R}^2$ to $\\mathbb{R}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e27d61b",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = np.linspace(-1,1,100)\n",
    "gridx = grid[:,np.newaxis]\n",
    "gridy = grid[np.newaxis,:]\n",
    "sq = - gridx**2 - gridy**2\n",
    "volcano = np.exp(sq) - 0.7 * np.exp(sq*4)\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "ax.plot_surface(gridx, gridy, volcano, cmap=matplotlib.cm.coolwarm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "330f351d",
   "metadata": {},
   "source": [
    "Now we have a function, we can build a filtered [cubical complex](https://gudhi.inria.fr/python/latest/cubical_complex_ref.html) from it and compute the persistence diagram of its **sub**levelsets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef58c852",
   "metadata": {},
   "outputs": [],
   "source": [
    "cplx = gd.CubicalComplex(top_dimensional_cells=volcano)\n",
    "diag = cplx.persistence()\n",
    "gd.plot_persistence_barcode(diag, legend=True)\n",
    "gd.plot_persistence_diagram(diag, legend=True)\n",
    "# red is dimension 0, blue is dimension 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9285415f",
   "metadata": {},
   "source": [
    "(The name top_dimensional_cells is because gudhi does the opposite of what we saw in class, it 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. The difference is not important here.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d98497f",
   "metadata": {},
   "source": [
    "Compare the 2 plots. Why do we only see 3 red points but 5 red bars? Why are more than 60 of the bars invisible? You can print `diag` to help, the format is a list of `(dimension, (birth, death))`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb9a4549",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7df4da5c",
   "metadata": {},
   "source": [
    "Now compute the persistence diagram of the **super**levelsets of this function (hint: there is no direct function to that, only sublevelsets)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec384536",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "be202de6",
   "metadata": {},
   "source": [
    "What happened to the point corresponding to the crater of the volcano between the sub- and super-levelsets?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eff023f",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9fc35f52",
   "metadata": {},
   "source": [
    "### 1d function\n",
    "In class, we looked at the function $f: t \\mapsto sin(t)+sin(2t)$\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": "code",
   "execution_count": null,
   "id": "32cf2ca5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fd65480",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c209c146",
   "metadata": {},
   "source": [
    "What happens if we consider a longer range instead of `[0, 2π]`?\n",
    "\n",
    "We will reuse this function later."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "114f66de",
   "metadata": {},
   "source": [
    "## Point sets\n",
    "### Torus\n",
    "As in the class, we first consider a set of points regularly spaced along a curve drawn on a torus. For simplicity, we embed this torus in $\\mathbb{R}^4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3258d48a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = np.linspace(0, 2*np.pi, 50, False)\n",
    "b = np.stack((np.cos(a),np.sin(a),np.cos(5*a),np.sin(5*a)),axis=-1)\n",
    "# Plot the points on the unwrapped torus\n",
    "plt.figure()\n",
    "plt.scatter(a, 5*a % (2*np.pi))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b188a6a",
   "metadata": {},
   "source": [
    "We now compute the persistence of the Čech filtration of these points. We actually use an α-complex for that. Notice that the data-structure used to represent a simplicial complex in Gudhi is called SimplexTree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e73caf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "cplx = gd.AlphaComplex(points=b).create_simplex_tree()\n",
    "p = cplx.persistence()\n",
    "# print only the most persistent features\n",
    "print([(dim,(birth,death)) for (dim,(birth,death)) in p if death - birth > .1])\n",
    "gd.plot_persistence_diagram(p, legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc0ba4a4",
   "metadata": {},
   "source": [
    "Can you recognize the features of a torus here? Is there anything extra?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8967c09",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b873a534",
   "metadata": {},
   "source": [
    "Now try doing the same computation, but instead of using an α-complex we will approximate the Čech complex with a Rips complex ([doc](https://gudhi.inria.fr/python/latest/rips_complex_ref.html)). What happens if you only change the name of the class? While the α-complex naturally has the ambient dimension, the Rips complex may be built up to an arbitrary dimension, so you need to specify a `max_dimension`. Can you still see the 3-sphere? What happens if you specify a larger dimension, say `max_dimension=5`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00eca1c0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "dc74ad04",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "75ce961b",
   "metadata": {},
   "source": [
    "Once the diagram has been computed with a call to [`persistence()`](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html#gudhi.SimplexTree.persistence) (or just [`compute_persistence()`](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html#gudhi.SimplexTree.compute_persistence) if you do not need the diagram in the form `persistence()` returns), you can get the points of the persistence diagram of dimension i as a convenient (n,2) numpy array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3876b61d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cplx.persistence_intervals_in_dimension(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c54072ab",
   "metadata": {},
   "source": [
    "### Time series\n",
    "Let us go back to the function $f$ define above. We saw in class that we can turn ([doc](https://gudhi.inria.fr/python/latest/point_cloud.html#time-delay-embedding)) it into a 2d point cloud with nice loops, so we try that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca6d04dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gudhi.point_cloud.timedelay import TimeDelayEmbedding\n",
    "f2 = TimeDelayEmbedding(dim=2)(f)\n",
    "plt.figure()\n",
    "plt.scatter(f2[:,0],f2[:,1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc2ecacb",
   "metadata": {},
   "source": [
    "Hmm, those loops are way too squished to see anything. Can you fix it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "801d2224",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "33f23129",
   "metadata": {},
   "source": [
    "Once the figure looks nice, compute or approximate the persistence diagram of the Čech filtration of this point set (here you have several choices). Dimension 1 seems the most relevant. Does the result have any connection with the diagram of the sublevelset computed at the beginning?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76aa72bc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0f2c2851",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0136687f",
   "metadata": {},
   "source": [
    "## Manually creating a filtration\n",
    "If you are not satisfied with the existing filtrations (AlphaComplex, RipsComplex, etc), you can also construct a simplicial complex by hand, specifying each simplex and its filtration value.\n",
    "\n",
    "Create an empty simplicial complex ([doc](https://gudhi.inria.fr/python/latest/simplex_tree_ref.html)) and `insert` a few simplices. You can see the list of simplices in your complex using `list(cplx.get_simplices())`. Notice that when you `insert` a simplex with filtration value `f`, the library helpfully ensures that the faces of this simplex are also present with a filtration value at most `f`. The function `assign_filtration` can be useful to change the filtration value of a simplex, but it does not provide this safety net."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74146547",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {
    "TorusTriangle.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAGSCAIAAABc3DAjAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd0AU174H8JndZanSeweVEoogKiCilGgs6DU0Yy+IYANNXnLvTW7ykpdy896NBlSwiw2liDX2jh0FiYjYKAqK0qvU3X1/zL2TDSDbZuZMOZ+/hnX3zDcb+O3szJnfQSUSCQJBnNfT05ORkdHZ2VlWVmZubp6YmAg6EaQSFJY2mtu8ebO1tfXKlStfvnwJOgubHT58OCwsTCgUIgjy3Xff2djYLFmyBHQoSHmwtDEDisL/U+Q6fPhweHg4tl1dXT1//vwLFy6AjQSpQiDPk3bv3q2vr8/n8w0NDQMCAsjOBEHUy8rKwkublpZWa2sr2DxcRkjB4cl8xrFjx4yMjGbOnDl9+vSSkpLnz58rtycIorOMjAx8+8aNGxMnTgQYhsuIKjiyS1tWVlZwcDC2PWrUqCNHjii3JwhihN7e3uzs7M8//xx0EI4iquDILm03b97U0dHBtg0MDC5evKjcniCIEdLS0r777js9PT3QQTiKqIIju7RJMzAwePv2rXJ7giD6u337dlBQkK2tLeggEIKoVnBklzYe74/naGlpvXv3Trk9QazX1NQEOoJKKioqzMzMhg8fjiBIeXk56DgcRVTBkV3auru78e3W1lZs4g8E9ZeUlPTll18ydJJKa2urpqamg4MDgiBisfjGjRugE3EUUQVHdmnT0tLCt1taWgwNDZXbE6Q0sViMIAj9S8a3336rp6e3ePHi3t5e0FkUNmHCBBcXFwMDA319faFQeOnSJdCJOIqogiO7tPn4+DQ2NmLb9fX1fn5+yu0JUs7WrVvj4uIQBImPj9+2bRvoODL89a9/nTBhQnR0dEdHB+gsiikoKGhsbGxsbGxqaurt7d21axfoRBxFVMGRPcf97NmzZWVly5cvRxDkp59+mjp1qpeXl3I7gzji6NGjycnJhw4dMjIyAp0FYhiiCo5ct++kp6dbWlo2Njb29vZGR0crsRuIa65cufLll19mZWVZW1uDzgIxDCEFB96ZCJHl4cOHcXFxu3btcnZ2Bp0F4hzF5rUNore398WLF9JXNyCSSCSSysrK9vZ20EFkcHd3379//8KFC2/dugU6i/Kqq6uZPqmFKerr62tqaogajZjSduLECUtLS3t7eyMjo507dxIyJjSg/Px8JycnW1tbAwODr776CnQcGRwcHI4dO/bFF1+cOXMGdBaFVVVV+fr6WlpaGhoazp07t6enB3Qi1nr37t1f/vIXY2NjMzOzoKCguro6AgaVqKylpWXIkCH4gAKB4Pnz56oPCw3I1dVV+n/fqVOnQCeSrbm52cvL69GjR6CDKAZvBIJZt24d6ESs9Y9//EP6rY6JiVF9TALOtd2+fdvf31/6EQ8PDzMzMxWHhfrr7e29cuWK9CMODg5Dhw4FFEdeLS0tpaWlXl5eKIqCzqKAa9eudXV14T+ampp6enoCzMNiBQUFDQ0N+I8ffPBBcXGximPK1a9tcFZWVn0eWbFixaxZs1QfGepDLBY7ODhItxJbvXr1okWLwCWSrampKSoq6tq1a+bm5qCzKCY0NPT+/fv4jxERET/++CPAPCy2cuXKgwcP4j8Sc1Vd9QM/iUSCTULBhIWFPXjwgJBhof6Sk5PxYx8XF5eff/5ZJBKBDvVeYrE4IiLi3LlzoIMo4+TJkwLBvz/79fX1s7OzQSdireLiYrzViqam5vXr11Ufk7DJH+fPn4+Li1u/fv2MGTOkb3CFCHf37t3s7GwPD4+oqKi2tjZdXV3a3ti7bt261tbWb7/9FnQQJV29ejUzM9PLyysyMrKnpweeZiHPmzdvvvnmm/r6+vXr19vZ2ak+IJHz2iZMmHD16lXpR3p7e/HPPYhr8vLyvv7661OnTvH5fNBZlNTV1dXe3g7vm6bGoUOHXr9+nZCQQMho5B5eZWRkPHv2jNRdQJimpiYCP6VU19jYmJCQkJaWxty6hiCIurp6/7pG4NwriDzklrZ58+bdvHmzpKSE1L1ACIJ0d3cnJSVhPUKAk0gksbGx33//vaWlJegsxMvPz+/z7QSiIdJPii1cuLCgoIAmf3KsceLEiT6PmJqazp07d8OGDXR4q9evX+/u7s7WlVOmTJkiEolyc3NBB4EGQ8X5/rlz58ILC8Sqr6/v/6Cpqeknn3wiPV8BiLy8vHPnzn399ddgYxCipqZmwN78ISEhvb29IpGI+kiQnOA5flYxNzcHO32sp6fn008/PXToEKNPsUl7323RISEhFCeBFELpwVRBQUFZWRmVe4QopqamdvbsWcbNzn0fFIWtcZiK0tI2cuTIixcvlpaWUrlTVpJ5B4JYLN6+fTuQP0ttbW3qd0oSgUCgoaEx+HOKiorgCpY0RPUpsNjY2IsXL7548YLi/XINj8ebMWNGcnIyPOhQhYGBgcwvnh4eHiiKnj17lppIkJwAnN1ftmyZiYkJ9fvlGjMzs+jo6LS0NNBB2C8kJEQgEBQWFoIOAv0BzGUE6UVrICXU19fLs+yApaXlkiVLKMgDhYaGgo4A/QnIORnd3d1VVVUAAzBX/3ltEBl6enpevXoFOgWkDJClTSgUnjx5Ei7TTY28vDx43k1RHR0deXl5Cr1ELBb36akHAQF4Jm1cXNy5c+dgdaOAra3thg0bYHVTCIqiit7dwePx1NTUjh8/TlIkSE7gbxKIi4urrq4GnYJhlFgjytzcfNasWdu2bYPVTX48Hk+JtysgIMDIyAhWN7BocTfC2LFjQUdgmD4d2+Vkbm4eFhbW0dEBL+PISSgUuru7K/HCgICAN2/eEJ4Hkh8tShtEmf7d3qFBqKmpubi4KPda1tySwVDgv5BK279/P7wgRQ3pBU0gUtXW1oKOwEX0Km3z5s07ffo0rG4yqX4Nrq6ubvPmzURkgWQoLS09evQo6BScQ6/ShiBITEzMb7/91tvbCzoIrVVUVKg4gpWVVVhYWEpKChFxWKu3t/f06dMqDuLn52dubg6rG8VoV9pQFI2Li4MrKlDAxsYmLCzs1q1boIPQmvT6mErDqltLS4vqQ0FyghWE0+zs7AhZPYitCGxq5OfnR8g4kJxod9QmraioCE55G1BERAToCJwA+7UxF61Lm7u7e3Z29uvXr0EHoZ0hQ4YQPmZ6ejrhYzIdj8ebP38+sWNWVlbm5OQQOybUH61LG4qiq1evxlYnBJ2F/QIDA1NTU0GnYD8bGxtHR8fMzEzQQViO1qUN+U91U1NTAx2EXt7XsF8Vtra2U6dO3bVrF+EjQ314e3sPGzYM3kVPKrqXNgRBUBSFrSv7OHDgABnD2tvbE/79i+kGXDxMdT4+PkFBQWSMDGEYUNpwIpGIkCvx0CDgAXIfsDUeQzGptPF4vH379sG7jqnx6NEj0BG44tKlS6AjsBCTShuKogkJCZmZmbC6WVpakr0LbW1teCcWNYyNjUk6w8BlTCptyH+uKjx48AB0EMAmTZpE9i7s7OymTp0Kq5sSrfEU5enp6e7uDiffEIt5dyPweDwK/rAhBEHs7OwmTZrU3Nysp6cHOgswyrXGU5Snpye8VkYs5pU2iEpDhw4FHYErLCwsQEdgFYZ9Ie0jMzOzrq4OdAoAioqKKN4jvN+IMtz8lSYcs0tbVFTU3r17OfirkJ+fT/EeX716xc3ZvNR3Rqmrq4MLY6uO2aWNx+MlJibu37+fjNn5kDRra+vAwEAOVrcnT55QvEcXF5exY8fC6qYixp9r4/P5a9asAZ2CE4YPHy6RSHJzc8ePHw86C8s5OzuLxeKXL1/a2tqCzsJUjC9t3ATqGrGTk5OTkxOQXXONq6sr6AjMxuwvpH08fvy4sbERdAoqUDBlF8LA1ngMxarSNnz48F27dsH7TKlx7Ngx0BGoQEZrPIW0trZy8BSn6lhV2rDzbmlpabC6UcDV1RWe6qbAkCFDQkJCtmzZAjoIw7CqtCEIwufzExMTWb++xsGDB0FHQJycnAICAlhf3ehw8d3e3v7DDz88cuQI6CBMwsLLCAKBwN7eHnQKctFkgWQnJyfWr5F+4MCBRYsWgU6BDBs2bNiwYaBTMAnbjtr66OzsBB2B5XR1dUFHgKABsLm0iUSiTZs2NTc3gw5CPODntvurrKwEHYErLl++DDoCA7C5tPH5/LVr127fvp19M0JoOCPh3bt3u3fvBp2CeDScZzN06FDYbEomNpc25D9XFeD6GhRwdnb28/NjX3WjYQctW1vbadOmweXHBsfCywh9qKmpffzxx6BTcIKLi4tIJKqrqzM2NgadheVsbW3Dw8PFYjGPx/KjE6XB94WRKioqQEcYmJubG6xr1DA3N4d1bRDcemuOHj3a2toKOgUB4FdsylDfGk9R9fX1sJtef9wqbVOnTt2yZQs7qhvNVVdXs6PZP/Wt8RTV09OzYcMGWN364FZpEwqFiYmJW7ZsocmUVxazsLDw8fHJzMwEHYT9zM3NZ82aBatbH9wqbQiCCIXCtWvXqqurgw6iEh8fH9ARZHNxcXF1dYWLbFLA3Nw8KiqK/t+dqcT+K6T9CQSM/6/28PAAHUEunp6eoCOoioaTPwZkaWlJwyl4AHHuqE1aRUVFe3s76BQQrcF6wVCcLm2WlpYpKSltbW2gg3AC/GZKDZFIlJqaKhaLQQcBjNOlTSgUrlmzJjU1lXHV7cSJE6AjKMzY2JgOvZhYj8/nR0ZGbtiwgePVjdOlDUEQoVCYkJDw8uVL0EEUU19fDzqCwjw9Pd3c3BhX3RgXGEEQU1PTOXPm7Nu3D3QQkBh/Ql11GhoaH3zwAegUnODp6cm4exUYOk/I1NR04cKFoFOAxPWjtj44fgxPAXhWHqIGLG1/EIvF69ev7+joAB1ENjr0fVURUzpN0bA1nqJyc3NFIhHoFFSDpe0PPB5v5cqVGzduZER1Y7oXL15kZGSATiEbDVvjKcrT0zMpKam3txd0EErB0vYnmpqaq1evhutrUMDLy8vFxQXeiUUBfX39mJiY5ORkTlU3WNr60tTUnDNnDugUMjQ0NLDghkEvL69hw4a9efMGdBD209fXX7JkSU9PD+gg1IGljZFOnz7NjrMnPj4+NF8Ti7at8RRlYGCgqakJOgV1YGkbzJkzZ+h57R9FURYctTEC+1rjtbS0cOGbKSxtgxk/fnxycjINV/xDUZRl81Tq6upycnJAp+AEiUSSlJTE+i+nsLQNRktLa9WqVRs2bKBbdRs+fDjLmkcbGxs7OjpmZWWBDsJ+enp6sbGxycnJ7K5urPrzIANW3dTU1EAH+ZNRo0bRLZLqvL29HRwczp07BzrInzCiNZ6i9PT0li1bdvfuXdBBSARvtJJNS0sLdASuGD16NN3OITKlNZ6idHV1x44dCzoFieBRmwLevn3b3d0NOgXLoSgKOgLEBrC0KUBHRycpKYkO10zv379Pt9N/hGP31yVa2bp1K/s+s2FpU4C2tvbKlSuTk5OBV7fy8nL2/S72IRaLjx49CjoFI1vjKWr+/Pn0nAmgCljaFKOtrR0fH19cXAw2Bo/Ho9s5KcL5+vpaWFgAr25MbI2nKC0trRUrVuzcuRN0ECLBywgK09XVHTlyJNgMHJmy6+vr++TJE9ApOAH7RgI6BZFgaWOkiRMnamhogE5BBWdnZ9ARIEaCX0iVJ5FINm7cCOSeFS0tLZZN2ZUJ1PlNFrTGU1ReXh4L+npx68+DWCiKLliwgAv3rNBBYWEhF87o04GHh8e6deuYvo4lLG0qwe5Z2b9/P8X77enp4cK5Nmm+vr4mJibHjh0DHYT9NDU1P//8cyau9CYNljZV6enpLV68mOKdXrhwoba2luKdAufn52dsbFxZWUnlTrlwhbQ/dXX1VatWAZ/kpApY2hiJI1dI+wsICLCxsaFyj5z9FqypqWlkZAQ6hfJgaSMSZetrsK+pEURbXV1dTDzvBksbkTw9PdevX0/BNVNjY2OhUEj2Xuisubn55MmToFNwAoqimzdvbm1tBR1EMbC0EUlfXz82NpaC1YN8fHwY/WVBdXp6egYGBhR8W4QT64RCYUJCwtatW5lV3WBpIxi2vgbrb/Ckg7Fjx+rr65Pd383f35/U8RkBq255eXmggygA3o1APENDQ9ARuCIwMBBOKqSGUCgMDQ0FnUIB8KiNRK2trSSd7H/58mVzczMZIzMO+7oNQ4SApY1Evb2969evJ+OaaXl5eWNjI+HDMlpJSQkZw966dYuMYRlt37599P9khaWNRAYGBjExMZs3byb82I2z89oG0dDQQMZVBdh6pL85c+bs3LmzoaEBdJDBwNJGLgMDg1mzZt25c4fYYWFp6y8gIEBXV/fMmTOgg7Afn89PSEjYu3cvnSdXwssIpDMxMTExMSF2TFdXV440NVLIhAkT8vPzQafgBIFAsGbNGtApBgNLGyMZGxuDjkBThK+tFxERQeyAEDXgF1JK7dy5k87H8GxC1Bf2IUOGEDIOiz18+JCG591gaaPU9OnTk5OTYXWjwJ07d+i2WjNbubq67tmzp66uDnSQP4GljVKmpqZz5szZsWOHiuPk5+fDK3eD8/PzU1NTO3v2rIrjwBtLZMKuKuzbt49W1Q2WNqqZmZnFxsaqOIhIJKKmxQijBQcHCwSCsrIyVQY5cOAAUXlYjM/nJyYm0upjAF5GAED1BdK5sFgfIZh1bxCj8Xg8S0tL0Cn+AI/aQCooKFCuQmloaMAbjCB6kkgkTU1NoFPA0gaUlZVVcnKyEtXN3d3dycmJjEhs1d7efunSJSVeSKsjEabYt29fdXU12AywtIFkZmY2e/Zs5aobpBBtbW0EQZS4Zjpp0iQS4rAZiqKrVq06dOgQ2OoGSxtgZmZm0dHRLS0toIOwX0hIiEAgOH/+POgg7IdVt8LCQoAZ4GUE8JT4ytPc3CwWiw0MDMjIw2IhISHwU4QaKIpOmTIFYAB41EYj8rcdf/nyZXl5Oalh2EpXV1eh5xcVFZGUBCIVLG00UldXl5qaKs8z4YpWqquqqpLnafB+e9UdOXLk9evXFO8UljYaMTc3nzFjRkpKisxnwnltqnvy5MmFCxdAp+CEmTNnHjlyRM7PEqLA0kYv1tbWM2bMyM3NHfxpZmZmVlZW1ERiq9DQUBRFYXWjAIqiK1asOHbsWFtbG2U7hZcRaMfGxkbmAukcX6mPKKGhoVevXpVIJIPcHxIUFERhItZCUXTlypVU7hGWNojTJkyYMPgT7O3tKQkCEQx+IaW1zMxM0BEgiGAvXrx48eIF2XuBpY3W/P39B7xm+urVq+LiYurzsNjdu3dlnuKECGFnZ3f27FkVO7LIBEsbrdna2k6dOnX79u19Hu/s7ISL9RFr9OjR7e3t/avbwYMHgeRht2XLll28eLG0tJS8XcDSRnf29vbz58/v8yBc0YoMU6ZMaWtr69Pjs6urC1QedouNjRUKheSN/8dlhNra2szMTOxSUW1tbUJCgqGhIXk7huTXf/EqWNpIMnXqVNAROETmTABpihaoP0rbmTNnVq5cib2yqKgoLi4uOztb2cwQKZ49ezZ8+HAEQRwcHBwcHEDH4QS47AsFmpqa9PX1B3+OogXqjy+k2dnZ+OyeoUOHAm+3BPUnFAo3b94MOgUndHZ23rp1C4GL9VHi2LFjJSUlgz9H4QIl+Y/W1lZ8e8eOHQcPHpTIraCgICAggMfjjRgx4sKFC/K/EFLUzZs33d3dtbW1hw4dunfvXtBx2Cw7O3vy5Mn6+voWFhY//PAD6Dhs1tvbO23atCFDhmhpaS1fvryjo6P/cxQtUIj0D3l5eZmZmVu3br106ZL8sTo6OqRv+hkyZMjr16/lfzmkkHHjxuFvNY/Hu3XrFuhErNVndZ49e/aATsRav/76q/Rb/de//nXApylUoP50N8Lo0aO9vLxu376dk5MzevRoHR0dGYd8CIIgyMOHD1+9eoX/2Nraev369aioKARBWlpa8vLysMetra1dXFyw7eLiYvx40t/fH+uA2tvbe+XKFexBQ0PDkSNHYtvl5eX4ReIRI0aYmJhg27m5udgKOhoaGvgf/Nu3b/EuNMOHD7ezs8O27927h/drx24exKLeuXOnf7xHjx7hjQrweCKR6PLly/3jVVRUPH/+vH+8a9euYRfXpOPV1NQ8ePBg8HghISE8Hg9BkLa2ttu3b2MPWllZubq6tra23rhxA3+rxWLxb7/95ufn1yeegYEBvoi6dDxPT09TU9M+8dTV1QMDA/vHGzZsGD4RPz8/H59rEhwczOfzB4yHbZeUlOC/D35+fthvkVgsxvt3S8d78eLFs2fP+se7fv16Z2dnn3i1tbW///57/3gFBQX4Er94vPb2duwbJYIglpaWH3zwAbb9+PFj/D5tX19f7FSadDx9ff1Ro0YhCHLmzBlEypEjRxYsWNAnnlAoHD9+fP94Q4cOxU+GSscLCgoSCAR94llYWLi5ufWPN2bMGKwFk0QiuXjxIvagnp7e6NGjse2XL18+ffoU2/bw8DAzM8O2b9y40dHR0SdeXV0d3htSOt79+/fr6+v7xHv37t3Nmzf7x3vy5EllZWWfeAiC4LfiSserrKzELze7u7ubm5tj2zdv3nz37h2CIGpqatjdIH3e6tOnT//8889IPwoVqL43WqmpqQUGBnZ0dPz0008//fTTIK/E6enp9XmkqKgIK23q6uqOjo7Yg9J9skxNTTU1NbFt/AIwj8fDn4z/K4IgBgYG+ONYlcHY2dlhC9Zh/zMwOjo6+JOlT0xaWlr2v54iFAoHjGdiYoJflMRXV0FRdMB4+vr6+ONaWlr447a2tv3jaWtr40+Wft+k4+EnFPrH09TUVFdXx/6oMPjNpO9796Tj9Xn3sPZw74sn/e5ZWFjgbS+xstsnnvS5dhMTE3V1dfw5+H8U/mTpC756enoDxrO1tcXiYXUKo6Wl9b54+I/4u6empjZgPGNjYzzV4PH09PTwP2PshYPHk373pLuESseT+e5Jx8PfRgRBFH33enp6EPnePXNzc/xXEY8n/e5JVxAjIyP8L2LAeNIP6urqDjiIjY1Nn3h9asgglxQUKFADHsu1traOGjVKzoNJsVj80Ucf4QOOGDGis7NTztdCikpMTMTfamNj419//RV0ItbasmUL/lYLhcIrV66ATsRaly9fxqsqgiBZWVmDP1+eAvXHcLNmzcK3W1pa5G9OjaLo0aNHk5KSzMzMfvzxx6tXr0pXbohY69ev37Nnz5w5c5YuXXr//v1p06bt2rULdCh2iouL++2338aPH5+QkJCfn4+daQEdip2CgoJu3749efLkgICAS5cuYd/5+lC4QOFFrrKyEt/etGnTgQMHFC2948eP7/PIjRs3FB2Ey2bPnq3Eqx4/fnzx4kXCw0CYBw8e4NsnTpy4fv06wDDslp2dja3uNiBFC9Qf51k0NDSys7MbGxu7uroEAsHs2bNVL8bGxsZpaWmLFy9WfSjWKygoOHjw4IEDBxR9obOzs7OzMxmRIARBPDw88O2wsLDTp0+LRCLpc1gQNRQtUH+UNmNj4wGPA1WBLQMMq5tMLS0tcK0D8kyfPv2XX34xNDRMT08PDg4eMWKE0kOBXaWJyxQtUKTfHu/k5BQcHAz/bgd3584dbMKBnNrb2/F5BtJOnz5NXCj2ePz4sa+vb0BAgJ2dnSp1DWIQKjp/2NvbS1/+gPp48OABPm9ITt3d3QM283N0dNy9ezcxsVjkhx9+aGpqevr06ccff6zoa9+34HxhYSE+ZxOiIUorDjaNEJLW1dX19u1bRVdZfl/nD2dnZ39/f1jd+qurqysqKhrwTRvc+1aZ8/LyqqqqgtWNtigtbeXl5enp6VTukf7OnTsXGhqq6KsGaWrk7Ow8ffp0lXOxyuXLl589e2ZhYfH9999j8+AJER4eXlVV9ejRI6IGhAhE6bIvH3zwAY/H27t3L37DCse9fPny3LlzL1++xB9JSUkZPnz4pEmTBn/hkCFDJk6c+L5/hetd9fHFF19g0+InTJiwefPmzz77jKiRw8PDlTgShChAZEfDCRMmXL16VebTioqKhEIhnK/QH+ENJt++fYvfVwghCFJbWzt79my49ig9HTp06PXr1wkJCYSMBuDsvoeHB6xr/WFFjdjS1tDQsGfPHgIHZKLm5mb8XdXX18dvBSdWT08PfvM5RAfwwiUt3Lt3D/uwSkhIuHfvnsznSyQSrOvJ4FxdXUeOHKnENGA2CQsLw0/2V1ZW4h1KiKWmplZaWnr37l0yBoeUALK0dXZ2Hjp0CGAA+hg1atTGjRslEsnGjRvlmeDW1dV15MgReUb28PBwd3d/+/atyhmZasWKFd7e3tj2vn371q5dq9DL5T/Ki4iIePnypTyfTBAFQK4er6GhMWzYsMzMTOkbXyF5KHRWztPTk9QwNDdjxoyNGzfy+XyBQDBy5Ei8m5icTpw4sWjRIjmfHBER8fjxY4UjQiQAWdoQBPHy8hKJREePHp05cybYJMwCV7SSn7a2NoGXRGXCG5pCYIE/1+bj4zNt2jTQKRiGx+MpcemzpqYmMzOTjDxQf/KcDIXIA760IVKdbCE5CQSCkJAQRV9lamrq4uKSkZFBRiS2Uvpq/vHjx+/fv09sGEh+tChtODixm2wjRoxwdnY+f/486CCM4e/vr9wLIyMjS0tL8/Pzic3DVp2dnQcPHsR67hOCsNJ27NgxfLkKpbW1tcl54Q9Smre39yB3MkAEioyMLC8vx9YBgAb32WefvXr16vjx40QNSExpKysr+7//+7+WlpZTp06pMs6YMWMsLS2PHj1KSCp2wxfugugsMjISnm+RKTs7u6ur67/+6786Ojr27dtHyJgElLaenp6lS5empqY6Ozv/+OOPZWVlqozm6+trb28P+7vJRMg3Heml/6AB4UvqQSQpLS1NSkrasGEDgiBRUVE7duwoLi5WfVgCSttnn322YMGCESNG8Pn8Xbt2LVy4UMXmCl5eXrC/GzW0tLSysrJAp6A1fKxW4B4AACAASURBVCVNFRUXF+NrvEK4rq6uJUuWbN++HVvoUk1Nbe/evbGxsS0tLSqOrGoFOXHiRH19PT6n0dnZeeXKlXFxcSoOi4FTt8jm7e09fPhwWN0o4Obm9vjx44KCAtBB6OXTTz+NiYnBF8BGEMTOzu7LL79ctGiRin/+KpW2Fy9e/PTTT9LLNSII8sknnxgYGGzfvl2VkTGFhYW//fab6uOwUlBQECHjeHt7KzpBH1JOdHT048eP4TQAHHaKrX+Ls7CwsA8++ODXX39VaXSlV9bq7u4ODg4uLCzEH8EX6+vu7v7www/z8vKUHhx38+bNY8eOqT4OBCmnrq6O2AG7urqIHZChnj9/Pnbs2Pb2dvwR6cX6RCJRWFiYKstaK3/Uhp9i6/9Pampq+/btW716dV1dnQpVF0EQxN/f38DAAN6XRw0CO9CyBuF9PYVCIbEDMlFXV9fixYt37NiBnWLrj8fjpaWllZeXK70LJUvb4cOHW1paBrlt2NzcPCkpqaurS8lcUgIDA+F9edR4+PAhnHlDDZFI9OzZM9ApgPn000+XLl3ap8FUe3u79I/GxsbyNyboT5nS9uLFi3/9618bN24c/Gl+fn5WVlZKpYJkO3jwIOFjjhkzxsLCAlY3aSTdCsrn8/Py8h4+fEjG4DR3/Pjxzs7OPqfYcnNzV65cOeA6bcpRprStXbt2x44dQ4YMISqEnDo7O2HrZxwhR8T9+fr6mpmZvXr1iozBmYi8Rp5z584tLCzkYHULCQnpc2BUW1v7+eeff//99zk5ObW1tYTsRZnSlpGRoei6mYTQ0NAQCARnz56lftec4u/vDw+3qTFv3rzOzk7QKaimo6MjfYpNLBYvXrz4l19+sbGxCQ8Pnz9/PiF3kipT2gY8D3r+/Pnq6mqV88gQFBQkFArft+otp1B/1AyRQZ6myuz2z3/+c9y4cYGBgQiC2NvbT548+euvv1Z9WGIm/VdVVX399df19fWXL18mZMBBBAcHw3lYCIJERESQvYuGhoYTJ06QvReaU3T1a6VJODlBPTc39/bt21988QX+yJo1a6qqqg4fPqziyASUtt7e3sWLF6ekpLi4uHz55Zeq9/+QycDAgOxdQAiCGBoaGhsbE9iMgYlkrglLlJycnJKSEmr2RRPYKba0tLQ+N1ampqb+61//UnHKFwGl7euvv545c6aPj49AINiyZcuCBQso6y9KwVdgjvP39zcyMoJnACgQGRl57949Qu4MZwT8FJuxsXGff9LR0dmzZ8/SpUtVmWipamk7ffp0WVnZypUrsR9HjBgxZ84cylrRP3r0iLNdFSsqKqjZUUBAQGhoKDX74rj58+fn5+dz5MLCzz//jJ9i68/JySkxMfGbb75RenyVSltVVdV33323bds26QeXLl3a0dFBzcq+oaGhPB6PmzNCrly5Qtm++Hw+ZfuiG4r74i1YsEBDQ4PKPQKRm5t769Yt6VNs/UVFRf3www9K70L50oafYtPT0+vzT5s2bdqxYwc1jeFDQ0P5fD6BfYehQXCw/yXs1UG4951i60+VKq98acNPsQ0YaP/+/StXrmxublZ6fPkFBwdz+bCCSq2trVy7ZoqiKJDGqE+fPiWqVRzdFBYWrlu3rv8pNmIpWdrOnDkjfYqtPzs7ux9//LGxsVHZYJAMA36okG3s2LGGhoacqm48Hg/ItAwnJ6c7d+6wsgPSxIkTx40bR/ZelCltVVVV3377bZ9TbP0FBwfb29srE0pZBQUF169fp3KPAHl4eADZb0BAgK2tLZBdAzF+/HgURYHsesGCBffu3ePajBCiKFPaEhISBjzFBtzIkSNbWlquXbsGOgjLDdjJiq1sbW0BtrNfsGABZXOGWUaZ/2f79u0D8m1IHlOnTm1qaoIfdNSAq/NQgIbHEIygTGnT1tYmPAeBpk+f3qcPFCvRYRrt7du36RCDCyQSyevXr0GnYBK4cBRT0eEXfezYsUKhkN29WE6dOtXR0QE6BYKi6KlTp54/fw46CGOwubR1dXXdvXsXdAqWCwoKEggEL1++BB2ELO3t7TT53r106dLLly+XlpaCDsIMbC5t6urqb9++hasIky00NJTF10xBTf4YUGxsLIs/RYjF5tKGIEhYWNjbt2/z8vJAByGeKm3jIfnp6OiAjvAnwcHBoCMwA8tLG4Ig4eHhJiYmoFNwQnNz88WLF0GnINhHH31Et+oGyYP9pQ1BEAcHB9AROEFPTw9FUc72YqHYkSNHVFnLjvU4UdpwbFpns76+HnSEAYSEhEgkEgqaLUMff/zxhQsXysrKQAehKW6VtpMnT7KmkQNtb+ScNGkSmzq8v379mrLWqoqKjY29cOFCn+U7IQy3SltUVFRZWVl+fj7oICzHppNTBQUFra2toFO817Jly2g+hR4UbpU2BEEiIyNfvnwJ+7tRg7JWwORBUZQ+kz8g+XGutCEI8vHHH7Ogv5uzszPoCLKVlJRcvXoVdAqVMKW0lZeXE7j0OgtwsbSxg7+/P+gIsk2ZMuXdu3eMrm4uLi7S6wHTloODw5kzZ+A1UxynS9v9+/c52BGbYlOmTGH0MbKjoyNTTmbFxcWdO3cOVjcMp0ubt7f3gwcPHj58CDoIy1HQUhXCxMXFqampgU5BC5wubQiCzJ07t7Cw8OnTp6CDKOzWrVugI0B0ZG1tDToCLXC9tCEIMm/ePEdHR9ApFMbENUFu3brFuA7v9+/fr62tBZ1CGdQsukRbsLQhCIIIBALQETjB39+fcR3eGxoaGLrm8cGDB1kw+UZpsLT9obu7+9mzZ6BTsBzW4Z1Bp7pp1dRIIfHx8WfOnOFsdYOl7Q9CofD27dtMWR5t+vTpoCMoafr06QxqWMCUeW0DiouL4+wcAFja/mT+/PlMWR7NyMgIdAROCAoKsrOzA51CSSiKMvcjUEWwtPW1YMECOvTC54K2trbbt2+DTgGxEyxtAxg5ciToCLLRth2F/HR0dOrq6hh3zZShTp48SYelgigDS9tg6HyS5cCBA6AjECAsLKypqYnO61e0tbV1dXWBTkGAqVOn5uTkcKe6wdI2mPT0dNjqj2xhYWF0nlf44MEDdqy0gqLoqlWrcnJyODLfDZa2wcybN+/ixYtweTSyWVhYgI7wXoy+QtoHVt04shw9LG0yxMbGXrt2jYb93SwtLUFHIF5DQwPoCH2hKEqTdUgJgaIo6AgUgaVNtkWLFtGwd8WkSZNARyDelStX6LYqtrm5OSsPc6qqqt68eQM6BYlgaYNoJDw8vLKyklbVzd7ens7fl5VmZWWVlZXF4qsKsLQpoLCwkB1nlOksPDy8traWNae3aAtF0dWrVx86dKi6uhp0FlLA0qYALy+vU6dO0eSmPBbfQDN16lTunBICCKtuLS0toIOQApY2xdDnlmO4Lhc1Kisrq6qqQKcgC4qijFhkQwmwtCksPj7ewMAAdApOyMvLA75ubH19PT1XsyYcC+5vkQZLmzJYecmMhsaMGQN83Vg2zWsb3JYtW2pqakCnIAwsbcrr6ekB+KsQFBQEatdUioyMLC0tff78OagA3Cltq1evPnjwIGtmhMDSpjyBQJCZmQnq8rm9vT2Q/VIvOjp66NChoPY+fPhwtp6N6gO7qnDlyhXQQYgBS5vysNtWWHz5nD4AXjDV1NRkxDqkhODxeJ988gnoFMSApU0l2Afd48ePQQfhhI6ODhZPeYGIBUubqlAUDQ4Opn6/OTk51O8ULE1NzZKSEoqvmYpEot7eXir3SBPnz5+vq6sDnUJ5sLQxVWtrK+gIAERHRz9+/JjKY7fnz5/fu3ePst3RR2ho6L59+5hb3WBpI1J6ejpD16xkkDlz5giFQsp2x50rpH3weLyEhIR9+/Y1NTWBzqIMWNqINGfOnOzsbGo+6IYMGULBXuiJykuWnC1tCILw+fyEhAQdHR3QQZQBSxuRUBSNj48/cOAABf3dIiIiyN4F/VHQ2ltLS0tXV5fsvdAWn89n6ALksLQRDDuMp2F/N1bKyckhe91YKysrd3d3UnfBCLW1tTRsFDoIWNogBpszZ05BQUFxcTHoIOxnaGi4e/duBt1OC0sbiR4+fNjY2EjS4HToPkIH8+bNKy4uZlOPb3ri8/mJiYkMqm6wtJHI1dV1165dJB3Gs+aGGNVFR0fzeGT9Jre0tMD+oxg+n79mzZoXL16ADiIXWNpIhP0qpKWlMeskBSStqakJLmmG4/P5jFiAHIGljWzYYfy7d+9AB+GEgoICwm96Y9mKVgSi+ZwYWNpIJxAIrK2tCR/Wx8eH8DGZbuTIkXl5ecReM+XxeDT/GwYlJSWFzrN5YWmjjkQiIfDwzcPDg6ih2GTBggV5eXkEfoU0NDSEb/WAli9fvnPnTtqebIGljToikSglJYWb935SadGiRTY2NkSNpqmpaWZmRtRobIKdbDly5AjoIAODpY06AoFgzZo1W7ZsgdWNbFTeZMplAoEgJiYGdIqBwdJGKTU1tcTExKtXr6o+1Llz51QfhN26urrKy8tBp4DAgKWNakKhMCwsTPVxWLzuN1HU1dUvXryo4roKDQ0NcAqhPHJzc5ubm0Gn+AMsbRCbLV269PLly6pcVRCLxe3t7QRGYquAgIDt27eTd/uNomBpA+nQoUNwyhvZYmNjVfl743JTI4Xw+fy1a9fu27ePJsvRw9IG0owZMzZt2tTW1qbEa+fMmUN4HrYaNWqU0q8VCAQaGhoEhmExPp+/fPlymrS9gaUNJKFQmJCQsHnzZiW678OLgEpQ4vhLT0/vww8/JCMMK6mpqWlra4NOgSAIwsgmc2yioaHx6aef0uSDjvXS09PHjRtH4BKuPT09mzdvVldXFwqFDQ0Nra2tX3zxBXcW9xtES0sLiqIAe0HDozbwlKtrTOktQyvz5s07d+5cWVkZUQPGx8d7e3vHxcUtXrz4s88+8/Lymjt3LlGDM5qmpibYKZywtNHIs2fPOjs75XzyiRMnSA3DVsuWLbt8+bL8ZwB6e3urqqre9685OTnGxsb4j+PGjbt+/bqqEVlBTU0N7AR1WNpoxNraesOGDfJXN0g5MTEx8vf77+zsvHPnzvv+dfLkyXPmzHn48CH2486dO9esWUNARFbAJqj//vvvQPYOSxuNaGpqrlq1ClY3Whm8qdGuXbvGjh07evToNWvW7N69OyAg4KuvvqIyHs0JhcJx48YB2TUsbfSipaW1cuXKyspKmc+kcsE6tvr9999lvtWDNzXS0tL68ssv16xZU1BQ8Le//Y1W0/E5DpY22tHW1h4+fLjMp/n7+1MQht1GjBhx8uTJwVeZUFNTG6Sp0a1bt/72t7999913V65c+fzzzyMiInbt2kV8UFbYvn07lfd1wNJGa7C/K9ni4+PPnDkzSHUTCASurq7v+9fVq1f/4x//EAqFPB7vs88+S05O/vbbb8nIyQILFy5MTU1VboK6EmBpo6+enp5169Z1d3eDDsJy8fHx+vr6yr22oqLC3Nwc/zEyMpKC1bUZCpugvnv3bmp2B0sbfampqa1YsSIpKWnANdJv3bpFfSS2Urq0zZ07d8+ePfiPp06dio+PJygUC6mrq69atYqafcHSRmva2trx8fHHjx/v/09PnjyhPg+79fb21tTU9HlQLBZfuHDhfS9Zv359V1fXDz/8sH379qSkpI6Ojn/84x8kx4TkAm+0ojtdXd2oqCjQKTiBz+dnZmaGh4dbWVlJP15dXT3ISz7//HPyo7FQXl6eu7s7eTelwaM2CPo3FEVXrVp1+PBh6TafsKkRSby8vBRte9PT03P8+PFjx47985//zMzMHPzJf5Q2hV4GAXHy5En8DqHp06eDDcNKWHWT/rIPSxtJhEJhYmJiamqq/DNCkpOTp02b9pe//OXvf/97SkrK4E/+43/bL7/8snbtWuxW7fHjx+fm5iqadcKECYR0/Yfep7m5efv27WvWrJH/PiFIdd3d3bCFFEk6Ojra29ux+3APHTr0+vXrhISEAZ956tQpTU3N4OBg7MfW1lYZbUUkEolEIjl58uSlS5ck/9HS0iJRRFVV1dKlS3V0dGbPnv348WOFXgsppLKyMigoaNSoUWFhYdevXwcdh816enr+93//19/fPzQ09NChQ6DjsNzu3bvd3NwcHR03bNggEon6P2HmzJltbW1VVVVlZWXyDIgo9zJpPT09np6eeK20srJqbGxUdBBIThEREfhbrampWVJSAjoRa0lf60RR9PTp06ATsVZ6err08da6dev6PycwMPD69evPnz+vrKzcsGFDT0/P4GP++3tNfX19YWGhubm5urr6xo0bly9fLudXnufPn+/du/fBgwf4I69evbp27ZpIJMKvoy9btgzbqKmpOXr0KLbt5uYWEBCAbV++fPnZs2fY9qxZs/T09BAE6erqwmcMWVtbT506FdsuKCi4d+8etj158mRbW1tse+/evdhd5bq6up988gn2YGlp6cWLF7Ftf39//I6ZY8eOvX37FtuOjY1FURRBkNraWny92A8++AC/rffKlStPnz7FtqOjo7E5UN3d3fjkQysrq2nTpmHb9+/fv3v3Lrb90Ucf2dnZYdv79u3r6OhAEGTIkCGzZ8/uH8/Pzw//hDh+/PibN2+w7aVLl/J4PARB6urqMjMzpVe07ejowP5D8HNDUVFRBgYGCIL09PSkpaUNHm/SpEl4U8b9+/djqzTo6OjgrcnLysrweQ++vr4jRozAtk+cOIFfNIyJicFOYtTX1+fk5GAPurq6BgYGYttXr17F40VGRhoaGiII0tvbi9+QZGlpiS/xVVhYmJeXh21PnDjRwcEB205PT8fOyGhra+MN0crLy8+fPz94vCVLlmC/yQ0NDYcOHcIedHFxGT9+PLadm5v7+PFjbDsiIsLIyAhBEJFItHPnzq1bt+JvtUQiycjIsLCwwLuASMc7cOAAdjpcS0tr3rx52IMVFRX4gopjxozx8vLCtn/77Tf8MsXixYvV1NQQBGlsbMzOzsYedHZ2njBhArZ97dq1kpISbDs8PBz74iYWi3fs2IE9aG5uPmPGDGz7wYMHt2/fxrY//PBDR0fHPvE0NTXnz5+PPfjixYuzZ89i26NHj/b29sa2T548+erVK2x70aJF2DfxpqamrKws7EEnJ6egoCBs+/r1648ePZIzXlFRET4ZMzQ0dOjQodj2wYMHk5OTESkZGRmffvop8mfPnz+3srLCfmNHjRq1ffv25cuXI4PAKpyFhUV5eTm2ffPmzdTUVPkrrnRdw5w4cUL+l0Py6+jowP4ScP/85z9Bh2KtYcOGSb/VixYtAp2ItfDyhxkzZkz/53h6euLbHR0dkyZNGnzMf18hNTExwT/Avb298WMrebi6uuIfmNhQ8F4TkmhoaMycORP/UVNTs7GxEb7bJMEPrhEEQVEU/yoAEU76re7/I0a65adQKOw/ubovrMKFhITg1U4kEnl5eSlUdF+/fh0XFzdkyJD58+c/ffr09OnT165dU2gESE4tLS1///vf/fz8/P39b9682dDQsHHjxgFPu0Iq6unp+eWXXwIDA93d3Y8cOQI6Dsvt37/fw8Nj2LBhqampYrG4/xMWL15cX1+PbTc2NkZHRw8+IKLcywY0fvx4fLurq0uJESAl1NTUFBUVgU7BZjdv3sQ2iouLa2pqwIZhsezs7OTk5Pf96/3799evX49tp6enX7lyZfDR/v2FNCEhAT9nf+rUqRUrVqhyeInApeQoZGJi4u7uDjoFm+Gt8VxdXQ8cOIBf4YGo5OXlZWlpmZycvGfPHh0dHfxKy/sI8Jc9efIkOTlZX1/fwMBA5svk9/jxYxcXF6JGY4EnT56Ymprq6elh1z0hBkFRFFs3Njo6WvrUD0SNWbNmyf/kP/66Zs2alZiYuHDhwj5XK1RUVlZ27do1Agdkuh9//NHQ0JDP56MoiqIoj8e7f/++EuMUFRX1eaS3tzclJUUC7woiE4qi8fHxcPEK+iP9wGHq1KktLS1wBTOcn59fc3MzdjpALBZnZGTgU4oUkp+f3+cRgUAQFRWVlJQEe/MSq09rPB6PZ21tDSoMJCcqvhNNmzZNXV2dgh0xgr+/v66uLradm5uLT6YlhKmp6SeffJKRkUHgmND7WuNJJBLK2mFDiqLodM/o0aOp2RH94cdozc3NxsbG2traxI5vYWGB304AkUosFm/btq2hoQF0EGgAVJ/JhmeCcEVFRW5ubkq/HL/ZBQKFz+cnJibu3r27vr4edBaoL6pL24ULF/CbBLlMJBINuOKB/PC7RwaRm5sLP0tUN0hrPKy63bx5k8o8kDyoLm0TJ06sqqqC1e3YsWP4zdXkcXR0lNmxD5IJu2f+ffh8PmwLSkMAplaFh4dXVlbW1tZSv2v6uHTpkvQibySxtraeOXMmrG4QB4GZNRoREWFiYgJk1zRRXV2t4oIXeAehwVlbW4eHh/f09KiyL0hOOTk5TU1NoFNACAKXfQGlrq5OxRFaW1vlfKaFhUWfVkiQQuRfFXjmzJk7d+6E10zpAHBpq66uVm4uPtMNfvqGJHASFtmwqwppaWkqXiOCVAe4tFlYWJSWlvafWM96hw8fVnEEGWteDKSurm7btm0q7hcanEAgWLt2LZyjDhz4L6SRkZFlZWUFBQWggzCM9CIJcrK3t588efLmzZvJyMNilpaWCj0fNj6gA1r8P4iKisLXEIBIZWtr+9FHH/Xv+Q4NYtKkScq9sLS0tKWlhdgwkJxoUdoQQOeeuMnR0VF6BTKIPLa2tlu3boXVDQi6lDYcvtAUNLiKigrQESAZ1NTU1qxZs3XrVvkvZ0NEoV1pu3DhwsOHD0GnYIArV66oOEJvb+/+/fuJyMJy/VvjyU9NTS0xMRFfPBCiDO1K29y5c3///ffi4mLQQdhPIBD4+vri64FC76PiFXyhUOjk5ERUGEhOtCttCILMnTsXX+QVItXw4cPHjh178OBB0EG4ore3F3QErqBjaUNUuCbFHT4+PoSM4+LiEh0dTchQ0ODEYnFSUlJ7ezvoIJxA09IGyeTh4UHUUHw+n6ihWImo1ng8Hi8hISElJQXeFkIBupe2c+fOPX36FHQKroDz3QYkT2s8OQmFwsTExJMnTxI1IPQ+dC9tkyZNun37dklJCeggnCAUCvHlaCGSqKurK7ToHKQcupc2BEEWLFhw8+ZNjvd36+/cuXOEj+ni4uLr6wurG8QCDChtCILExMRwvL9bf69fvyZjWBcXlzFjxsD+btLkbI2nhJMnT757946kwTmOGaUNopKrqyvs7yaNvHsJJk6cuGnTJnhVgQwMK21v3rwpKysDnYIr4FLNZBMKhQkJCZs3b4bL0ROOYaXN3Nz8woULsLohCELBYqNPnjzJzMwkey/0R2rzNQ0NjYSEBKFQSN4uuIlhpQ1BkGXLlsHqhiAIBX8Mrq6uLi4ucC362bNnkzq+uro6bPFGOEa+ocuWLYPLa1JjxIgRTk5Ov//+O+ggnPDq1Sv4zZQojCxtCIIMHToUdATAGhoaqKnvI0eOHDFiBAU7ggwMDDZs2ACrGyGYWtpwnP09OHv2LJyiQQ3KWuNpaWmtWrUKVjdCML60paWlVVVVgU4BAIqiFH8rF4lEJ06coHKPNKF6azz5aWlprVy58smTJ5Ttka0YX9qWL19+5swZDjZBor608fl8a2vr7OxsKnfKQdra2vAMgOoYX9oQBImJibl58yboFFQbOnQo9ZfVvL29bW1tjx49SvF+IUhRbChtKIpGRUWBTkG1UaNGAZkM5evrGxoaSv1+ASKqNZ6iJBJJUlJSd3c3kL0zHRtKG0QxJVZ3ZjQCW+MpBEXRZcuWJSUlwbXolcC20nbx4sU3b96ATsEVlZWVoCOwnJaW1ooVK+CsaSWwrbSFhIRkZWVxYQGhwsLCjo4OsBmqqqqOHz8ONgPr6ejoLFy4EHQK5mFbaUNRdPXq1VlZWQ0NDaCzkKuiogL49xR/f38TExPWVzcyWuNBZGNbaUP+U910dXVBByEXj8ejw91mWHVj96luklrjKeHy5ctwnracWFjaEATh8XgCgQB0CnJRP6/tffz9/WHjCmr4+PgkJSXB6iYPdpY2XG1tbV1dHegUpBg/fjzrj0yhPnR1dePi4mB1kwfLS5uRkdHevXtZWd309PTodmRaXFx85swZ0CmIR0FrPPlh1U0kEoEOQncsL208Hi8xMXHPnj319fWgs7Cfm5ubhobG2bNnQQchGN2+buvq6mpoaIBOQXcsL20IgvD5/DVr1rx8+RJ0EE4ICgoSCoWFhYWgg3BCY2Njb28v6BQ0xf7ShiAIn8/39vYGnYJgZ8+epefk5ODgYC8vL9ApiETbQ34ej5eUlASr24A4Udqk0eSqouroc4WU9WjbyklPTy82NhZWtwFxrrRt2rSpqakJdAoC0L+0iUSia9eugU7Bcnp6ejExMfn5+aCD0A7nStuKFSt27NjBgnsVjI2Nab5aKJ/P7+npOX/+POggLGdgYODr6ws6Be1wrrRhVxVY0HHM29vbxMQEdAoZQkJCeDwe06ubs7Mz6AiQwjhX2hAEEQgES5YsAZ2CK0JDQ11cXECnUIm/vz/oCPLaunUrnPKG4WJpgyhmY2MDOgJXzJo169dff4VXFRBY2nJzc1tbW0GnUEZlZWVzczPoFIppaWkBHYHl9PX1Fy9evGvXLtBBwON6afP399+yZQsTq1t5eTnjLobcuHEjNzcXdAqF3bp1C3QEBRgZGS1btgx0CvC4XtrU1NQSExOZWN1QFBWLxaBTKGbKlCnt7e2Mq25w6Twm4nppQxBEKBQmJCRQvzqUimjSr01RWHUD3kSTI+7cucO4zz+iMOzvmSTq6ura2tqgUyjGxcXFzMwMdAplTJkyRV1dHXQKTnB0dExOTuZmdYOl7U9aWlrevXsHOoVcjIyMuLayFCjTp08HHUFJJiYm8+bN42Z1g6XtTzQ092LDkQAADGpJREFUNDZt2tTe3g46CCeUlJQw4gy9kZER6AjKMzExmTt3bltbG+ggVIOl7U+EQmFiYuL27duBLxbFBa6urvX19Tdu3AAdhOVMTU052JAZlra+1NXV4+Li6H+/cVFR0dOnT0GnUFVYWFhtbe3vv/8OOggndHV1MfHSk3JgaRuApqbmuHHjQKeQoaurix0N8mfOnAlqeXY57d69G3QEYjQ3N2/YsIEj1Q2WNqZi4ry292HczBuGMjU1nT17dnJyMheqG/yVkmHbtm30nISlpqbGsoogEomKiopAp2A5U1PTWbNmXb9+HXQQ0rHqb4MM2LXzzs5O0EH68vT0dHNzA52CSHw+/+nTp/fu3QMdpC9LS0vQEYhkYWERGBgIOgXpYGmTQUtLC+teCToIJ0RERFRUVNCtuk2aNAl0BEhhsLTJpqOjs2rVKtApuCIyMhLeq0CZ/fv3s/W8GyxtTNXc3NzY2Ag6BSlofsGUTT788MNt27axsrrB0qaYvLw8mky5qKqqKi0tBZ2CXDRpGMviixvm5uZhYWHbtm0DHYR4sLQpxs3NLTU1lQ7VjU2TP97n6NGjDx48AJ0Cof/8bVVYWVnFxMSATkE8WNoUo62tvXjx4qSkJODVjaFNjRQSERHx8OFDFh800YRAIAAdgXiwtClMV1d32bJlwHthm5qacmHNgTlz5hQXF9NzaiH7FBcXg45AGFjalKGnpwe8G4ShoSHL5lu9zyeffAL2mmlQUBDAvVNJV1c3JSUFdApiwNKmku7ubpqc6obIY29vDzoCRWxsbGbMmJGamgo6CAFgaVPJu3fvkpKSYHWjRklJSUlJCegULGdjYzNt2rT6+nrQQVQFS5tK9PX1Y2JigKz8WF1dzbXz666urvfu3Xv06BHoICxnZ2cH/HyL6mBpUxW28iP1qzR1dHSwdcruIObPn5+Xl0dxo7qcnBwqd0cfjJ5dBEsbAYyMjEJCQijeKRcmfwxo0aJFFF8/YdxCjkSpqalh7mxeWNqYCkVRbpY2BEF0dHRAR+AEc3PzyZMnb968GXQQZcDSRrADBw5QcxhvY2ND/1bApBKJRFVVVRTsiMu369va2n700UeXL18GHURhsLQRbOLEidSsjcbj8Vg5iVx+fD7/9OnTFNxIO3v2bLJ3QWeOjo7BwcGgUygMljaCYWujsWbeI83FxsZeunSJ9W0CICXA0kY8U1PT+Ph4svciEongioIIgsTGxr569Qp0Cq44duwY6AjygqWNFGpqamTvoqam5tq1a2TvhRHGjx9P6vgVFRWkjs8gbm5uO3fuBJ1CLrC0kaukpISk65hcvkJKsStXroCOQBfDhg0bP378rl27QAeRDZY2chkaGpK0NhoX+rUpKjs7+8WLF6BTsNzw4cMjIyNBp5ANljZymZmZkbTyo5aWlrW1NbFjMl1UVNTp06fLy8tBB2E5XV1d0BFkg6WNdGZmZlFRUW/fviV22CFDhsA1BPqLj48/e/YssYsr+vj4EDgam1RWVoKO8F6wtFHBysrK3NwcdAquiI+P19DQIHBA+BHyPh0dHWlpaaBTDAyWNkrBE/8Qmzg5OY0dO5ae1Q2WNkrV19cT1eevo6OjsLCQkKHY6smTJ69fvwadguWcnZ39/PxqampAB+kLljZKGRsbT58+nZB7Fbq7u+F8q8E5OTnl5OSoXt3OnTtHSB62cnV1NTU1BZ2iL1jaqGZjYxMWFqb6Xwuc1yYTiqKrVq1SvbrBQz8mgqUNADs7u0mTJqk4CCxt8sCqG+gUXFFTU7N//37QKf4Nljam0tbWZmI/BuqhKMqRpb+AMzU1HTVq1J49e0AHQRBY2oD77bfflHshj8czMDAgNgy7SSQS5VaPnTNnDuFh2MrFxcXHx+f8+fOgg8DSBpqHh8fu3btBp+CKtLS0N2/eKPoqoVBIRhi2cnd3nzhxIugUsLSBZmdnN27cOKZ0U2A0FEUTEhIyMzOVqG4Q48DSBt6wYcPCw8MVfVV3d/fBgwfJyMNiKIquXr363r17Cr2Ks8u+qA5g53FY2mhBibNm8Aqpcng8XlhYmEIv4exifaozMjIC9QEMSxu9VFdXy/lMWNog+vP09HRzcwNS3WBpo5eWlhY5ryrweDwzMzOS47DckSNHGhoaQKdgOU9PT7LbIA8IljZ6cXZ29vf3l6e68Xi8Dz/8kPxEbDZjxozdu3fX19cP/jRnZ2dq8rCVlZUV9TuFpY12nJ2dR48eLf83U0hpfD4/MTFx9+7dg/d38/f3pywSuzU3N1O2L1ja6MjNzc3CwgJ0CmWIxeIVK1Y0NjaCDiIvPp+/du1aYvu7Qe9TXl6ekZFBzb5gaWOwoqIi0BH64vF40dHR06ZNo2Zdd0LwePCvgCJeXl4uLi6ZmZkU7Av+T6W1urq6Qa4u5efnUxlGTkFBQdu2bZs1a9aTJ09AZ1FMWVnZgN+Ybt26RX0YtvLy8ho2bBgFzVRgaaM1Y2NjNze3AwcOgA6iGHd39/379y9atIhZRcHOzm779u39v00zrkbTnI+PDwUNC2BpoztsZhDjuiE6ODgcPXr0r3/965kzZ0BnkRd23m3v3r3K3UUP0QosbQwwYsSIAfu7BQUFUZ5FAWZmZqdOnUpJSWHQDWF8Pn/FihXwJlNq1NfXHz58mKTBYWljMHt7e9ARZNDR0cnJyTl+/Pgvv/wCOou81NTUnJycpB+ZPn06qDDsZmRk5ODgkJWVRcbg/G+//Vb1UTo6OjZu3Hj8+PGOjo4RI0bAS+nkyc3NvXTp0q+//nr79m0XFxc9PT3QiWTg8/kff/xxenr6+fPnJ0+ejKIo6EQKOHr06C+//HLjxg1ra2t47wcZLCws2tvbHz58WF9fv27duoKCAnt7e0dHRwKGlqhMLBZLt2fy8PDo6OhQfVhoQLNmzcLfahMTk9evX4NOJK///u//jo2NBZ1CAdKr8wiFwnv37oFOxFoXL16UnoKTkZGh+pgE3GL99OnTPneipKSk+Pr6qjgs1J9IJAoMDOzu7sYf+eyzz2bPng0wkvyKi4tTU1MJWc2LGtHR0WVlZfiP4eHhX375JcA8LPa3v/3twoUL+I+BgYG5ubkqjilQ8fUIgvS/nHT+/PmXL1+qPjLUR3d3d09Pj/Qjd+7cEQgI+J9Itq6urqysrI8//jg7Oxt0FnnV1dVJ/1hSUsKg8Mwi/RGCEHU/luoHfh0dHTY2NviAurq61dXVqg8LDUi6iQKPx7tz5w7oRLKJxeKIiIizZ8+CDqKYuLg46b+UvXv3gk7EWsnJydJv9d///nfVxySgtEkkksLCwgkTJujo6Pj4+Fy+fJmQMaEBVVZW/uUvf9HV1XV2dk5PTwcdRy7r16//5ptvQKdQWEtLS0xMjKGhoY2Nzc8//ww6DpuJRKKvvvrKwsLCxMRk9erVhJysh+0MIXLdvXv3q6++On36NJ/PB50F4hBY2iASNTU1TZ48+fDhw3AlUIhicp2B3r17t76+Pp/PNzQ0DAgIIDsTxA4SiWTp0qX/8z//A+sapBBCCo7suxGOHTtmZGQ0c+bM6dOnl5SUPH/+XLk9QVyTlJTk5uY24C1iEPQ+RBUc2aUtKysrODgY2x41atSRI0eU2xOkis2bN584ccLW1hZ0EHndvXv39OnT33zzDeggEMMQVXBkn2tzcHAoLy/Htl+8eBEXF8egXg4sw5RVrBh9iq2npycjI6Ozs7OsrMzc3DwxMRF0Im4hquAoNtvTwMDg7du3SuwG4pS2trZ169Yxsa4hCHLixIlZs2YJhUIEQb777rtdu3YtWbIEdCiOUqXgyP5CKn1vl5aW1rt375TbE8Qd1tbWjL7chNU1BEGWLVvGuD6gTEdUwZFd2qTvWGxtbcX/r0MQK0n32NHS0mptbQUYhoOIKjiyS5uWlha+3dLSYmhoqNyeIIgRpJdcunHjhnRXG4gCRBUc2aXNx8cH7xZfX1/v5+en3J4giFl6e3uzs7M///xz0EG4haiCI7u0LVy4EP8cO3PmDFNa6ECQitLS0r777jv6N/tkGaIKjlyTCdLT0y0tLRsbG3t7e6Ojo5XbE6Q6pkz+YIHbt28bGRkNHz4cdBAuIqTgwD8VxhCLxXw+XywWM6sHNxNVVFRIJBIHBwcEQcrLy7ENiFngsi/MsHXrVqx9WHx8/LZt20DHYbPW1lZNTU2snInF4hs3boBOBCkDHrVB0J+MHDkSmw0vkUja2toWLFiwa9cu0KEghcHSBkEQC/0/mQEYX7Ub0wUAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "id": "9fb8f204",
   "metadata": {},
   "source": [
    "The following code constructs a small torus.\n",
    "![TorusTriangle.png](attachment:TorusTriangle.png)\n",
    "In this example, we do not care about filtration values (everything is inserted at time 0), only the topology of the full complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4191b0d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "cplx = gd.SimplexTree()\n",
    "for i in range(3):\n",
    "    ii = (i + 1) % 3\n",
    "    for j in range(3):\n",
    "        jj = (j + 1) % 3\n",
    "        cplx.insert([3*i+j, 3*ii+jj, 3*ii+j])\n",
    "        cplx.insert([3*i+j, 3*ii+jj, 3*i+jj])\n",
    "print(list(cplx.get_simplices()))\n",
    "# Since homology in dim p uses simplices of dim p and p+1, persistence is by default computed up to dim-1.\n",
    "cplx.compute_persistence(persistence_dim_max=True)\n",
    "print(\"The number of holes of dimension [0, 1, 2] are\", cplx.betti_numbers())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baf37d3d",
   "metadata": {},
   "source": [
    "By gluing some \"spheres\" and \"circles\" (the boundary of a simplex for instance is a *topological* sphere), can you build a simplicial complex that has the same homology as the torus but looks nothing like a torus?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9afedb9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "fe3b66a2",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "67998d8b",
   "metadata": {},
   "source": [
    "## Rips filtration and higher dimensions\n",
    "In class, we saw that the Rips complex is not embedded, it is defined as an abstract simplicial complex. Here we will see that the higher dimensional simplices that appear can actually generate some topology."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46d5b8c1",
   "metadata": {},
   "source": [
    "Generate 20 points evenly spaced on a circle in ℝ². Build the Rips filtration on those points up to dimension 8, and plot its persistence diagram. What do you notice? You can try again with 21, 22, 23 or 24 points. You probably shouldn't try to increase the dimension or the number of points too much as it will quickly fill the memory on your computer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79d35023",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5668a9d3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8ceb5e1",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cafa6d21",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f6458bc6",
   "metadata": {},
   "source": [
    "## Distance and stability\n",
    "### Point sets\n",
    "Let us consider again the point set on a curve on a torus in $\\mathbb{R}^4$, and compute the persistence diagram of dimension 1 of its Rips filtration. Now perturb each point randomly by a small noise, and compute the persistence diagram of dimension 1 of these new points. Compute the [bottleneck distance](https://gudhi.inria.fr/python/latest/bottleneck_distance_user.html#gudhi.bottleneck_distance) between these diagrams. Retry it a few times, maybe also with dimension 0 or 2. Can you confirm the stability result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0af97d3d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1457e041",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3c2acb28",
   "metadata": {},
   "source": [
    "Now do the same experiment with the alpha-complex. Note that there is an even worse trap than for the Rips."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b9cbe72",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "a80ed5de",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "78b4cf0b",
   "metadata": {},
   "source": [
    "In class, we said that Rips and Čech filtrations have close persistence diagrams in log-scale. Can you illustrate that on this dataset?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38950638",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "85328da3",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7168da8b",
   "metadata": {},
   "source": [
    "### Functions\n",
    "Similarly, illustrate the stability property on one of the functions seen above (either the volcano, or the curve)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}