Tuto-GUDHI-Expected-persistence-diagrams.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Expected persistence diagrams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Theo Lacombe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A topological descriptor for a random process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial illustrates the notion of \"Expected Persistence Diagram\" introduced in a work of Chazal and Divol https://arxiv.org/pdf/1802.10457.pdf. In a nutshell, the expected persistence diagram encodes the topology of a *random* process as a *deterministic* measure. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that given an object $X$, say a point cloud embedded in the Euclidean space $\\mathbb{R}^d$, one can compute its persistence diagram $\\mathrm{Dgm}(X)$ which is a point cloud supported on a half-plane $\\Omega \\subset \\mathbb{R}^2$ (see this tutorial https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-persistence-diagrams.ipynb for an introduction to persistence diagrams).\n",
    "\n",
    "Now, consider that our point cloud is random. Here, we will consider that it consists of $n$ points sampled _i.i.d_ with respect to some law $\\xi$, thus denoted as $\\mathbb{X}_n = \\{X_1 \\cdots X_n\\}$, with $X_i \\sim \\xi$. To a *random* realization of $\\mathbb{X}_n$ corresponds a *random* persistence diagram. What can be said about this object? Can it be turn into a deterministic topological descriptor?\n",
    "\n",
    "As often, the simplest way to derive from a random object a deterministic summary is to consider its average behavior. Namely, for each compact $K$ included in the half-plane $\\Omega$, one can ask: _how many points of my (random) diagrams belong to $K$, on average?_. It leads to the following definition of **expected persistence diagram** of $\\mathbb{X}_n$, denoted by $\\mathbf{ED}(\\mathbb{X}_n)$:\n",
    "\n",
    "$$\\forall K \\subset \\Omega, \\mathbf{ED}(\\mathbb{X}_n)[K] := \\mathbb{E}({\\mathrm{Dgm}(\\mathbb{X}_n)}[K]),$$\n",
    "\n",
    "where $\\mathrm{Dgm}(\\mathbb{X}_n)[K]$ is, by definition, the (random) number of points in $\\mathrm{Dgm}(\\mathbb{X}_n)$ that belong to $K$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chazal and Divol proved (under mild assumptions) that $\\mathbb{ED}(\\mathbb{X}_n)$ is a deterministic measure which has a density with respect to the Lebesgue measure on $\\Omega$. Later, Divol and Lacombe (see https://arxiv.org/pdf/1901.03048.pdf) showed that $\\mathbf{ED}(\\mathbb{X}_n)$ is stable with respect to $\\xi$, which states that point clouds sampled from similar random processes must have similar expected persistence diagrams."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A numerical illustration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Theoretically computing $\\mathbf{ED}(\\mathbb{X}_n)$ cannot be done. In practice, one is more likely to have access to some observations $(\\mathbb{X}_n^{(k)})_k \\in \\mathbb{R}^{d \\times n \\times n_\\mathrm{obs}}$, $k = 1 \\dots n_\\mathrm{obs}$. As usual, the empirical (arithmetic) mean converges to the expectation, that is\n",
    "$$\\lim_{n_\\mathrm{obs} \\to \\infty} \\frac{1}{n_\\mathrm{obs}} \\sum_{k=1}^{n_\\mathrm{obs}} \\mathrm{Dgm}(\\mathbb{X}_n^{(k)}) = \\mathbf{ED}(\\mathbb{X}_n).$$\n",
    "\n",
    "Here, the sum of diagrams must be understood as a sum of measures. Briefly, it consists in taking the union of points in each diagram $\\mathrm{Dgm}(\\mathbb{X}_n^{(k)})$ (counted with multiplicity), and dividing the multiplicities (viewed as a weight) by $n_\\mathrm{obs}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Experiment overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will showcase expected diagrams in the following context:\n",
    "- Consider a point cloud $X$ with $n_\\mathrm{tot}$ points, which splits in $(1-p)n_\\mathrm{tot}$ points sampled on a nice shape $\\mathcal{M}$ (here, a 3D torus or a 1D circle), and $p n_\\mathrm{tot}$ outliers, where $p \\in (0,1)$ is the proportion of outliers. Let us call $X_\\mathrm{true}$ and $X_\\mathrm{out}$ these two point clouds.\n",
    "- Computing directly the (alpha complex) diagram of the whole point cloud $X$ won't reflect the topology of $\\mathcal{M}$ in general, due to the presence of outliers. However, if $p$ is small enough, one can expect that by sampling $n < n_\\mathrm{tot}$ points on $X$, we will avoid sampling outliers with high probability. We will thus recover a diagram close to $\\mathrm{Dgm}(X_\\mathrm{true})$ \"most of the time\", provided we sampled enough points.\n",
    "- Repeating this process, and then averaging, we obtain (an estimator of) $\\mathbf{ED}(\\mathbb{X}_n)$.\n",
    "- As the random processes $\\mathbb{X}_n$ and $(\\mathbb{X}_\\mathrm{true})_n$ (sampling $n$ points on $X$ and on $X_\\mathrm{true}$ respectively) are similar (in Wasserstein sense), so should be their expected diagrams. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: %matplotlib notebook enables iteractive 3D plot.\n",
    "#%matplotlib notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import gudhi as gd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# A collection of utils functions for this notebook\n",
    "import utils.utils_epd as uepd\n",
    "\n",
    "# Comment this if you don't want to use LaTeX in matplotlib rendering.\n",
    "from matplotlib import rc\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical illustration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of point on the shape: 9990\n",
      "Number of points corresponding to noise: 10\n",
      "Proportion of noise: 0.001\n",
      "Total number of points: 10000\n"
     ]
    }
   ],
   "source": [
    "nb_tot = 10000\n",
    "prop_noise = 0.001\n",
    "\n",
    "nb_points = int((1 - prop_noise) * nb_tot)\n",
    "nb_noise = nb_tot - nb_points\n",
    "\n",
    "print(\"Number of point on the shape:\", nb_points)\n",
    "print(\"Number of points corresponding to noise:\", nb_noise)\n",
    "print(\"Proportion of noise:\", prop_noise)\n",
    "print(\"Total number of points:\", nb_tot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X = sample_torus(nb_points, r1 = 3, r2 = 1)\n",
    "X = uepd.sample_circle(nb_points)\n",
    "X_noise = uepd.sample_noise(nb_noise, \n",
    "                       ndim = X.shape[1], \n",
    "                       scale = np.std(X),\n",
    "                       type=\"uniform\")\n",
    "X_tot = uepd.add_noise(X_noise, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD3CAYAAADhRcqHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeI0lEQVR4nO2d3W/cxrnGn5Fsyfqy1/JX6jROIhdob4oisoI4aIG2gHPZu6T5D9SrRnFaNE2bDzUpkDpA4ShAgTbnLwhw7tsi6UVxTpu6kXOdokdO0yROE9vatbTS6sPaORdDerkrcjlDcsgh9/kBC63ImeXw6+E7My/fV0gpQQghYQwV3QBCiLtQIAghkVAgCCGRUCAIIZFQIAghkVAgCCGRHCi6AXEcP35cPvDAA0U3g5DKcvXq1ZtSyhNh65wXiAceeADLy8tFN4OQyiKE+ChqHbsYhJBIKBCEkEgoEISQSDIZgxBCzEop349Y9ziABoBZKeVrUctIudjeBnZ3gWYTOHQI2NxUy5tN4OBBYGgIGB0FxseBvb3OsvHxYttNzEgtEEKICwAuATgXsm4WAKSU7wghZvz/e5dFiQvJn3odOHAAGBkBNjaA1VXg8GFgZ0cJQrsNfOUr6bZx7RowPKxEZnxcCcfYGHDnjtr26Kj6nxRPaoHwbvTViNVPAnjb+34NwAUAx0KWUSByptUCtrbU97099f2++/LZ9syMXrm1NSVMo6Pq75EjSlhIftie5qwBCIrHsYhlxDKNBiCEekofP150a/Q4fHj/slOngH/+U3VZpKSlYRsn/SCEEPMA5gHgzJkzBbemnLTbwK1bwMSE+lSFzz/fLxxffKG6RFNTqrtCssP24WwAmPa+1wDciljWhZTyTSnlnJRy7sSJUAcvEsL6OnDzJvDII8oUP3myWuIQxcmTQK2m9nluTh2D3d2iW1UNrFgQQoialLIB4C0Ac97iGQDveN/DlpEE7O4qYTjGjhoA4OpVwH+mDA2pQVZaFslJfdi8Kcs576/PnwDAn53wZjoaUsr3w5albcMgcvOmuuhHRigOUbTbHcvioYfUDA0xQ7gek3Jubk7yXQzF9raaenRVEIQAPvhADSDqzlQUwdqaml7ljIhCCHFVSjkXts7JQUrSzfa2ckYqis8/7/aDmJxUU6MTE2omwfdjOHiwU0dKVafZVMIxPKymVnd21EyK7wfxve8B//hHvvvjD3KePg3861/d7SbdUCAc5vZtZSLnyeqquvn9KdFTp5L/1sgIMD3d+T9s2vKDD7pnXFotdcMeOZJ8u7pcv67aODSkPEFHR+1vs2xQIBwkT2G4eVN5Lx48WJwb9NBQZ2DRb4OUqivQbitrI0xcsqLd7lhorVax1pprUCAcYmvLvuPPxoayDA4dUk9PlwmKgt9lqdeVoJw8aWeb/vHf2WHXA6BAOMHennqK2+BLX+r08aem7GwjL0ZGOl0eKVV3SEo7nqEjI51u1iBPkQ7wrrvBiRN2xKHRUBf39etKGMouDmFMT6sZHSnV/t6+ne3vS6m6N488ku3vlglaEAXRamXf5282lTVis7/uKv6gpm9ZZDkV/Pe/K2tid9eepecqA7a7xeNPC2ZJvZ7/bIfLTE+r49xsqgHIrGZE/IHcjY1sfq8MsIuRIz/9abbisL2tbgSKQziTk8qa8rsgWbC5qawJP0BO1aEFkQNZD0JubDAykylHjiihyMrpzH8JznFH5NTQgrDMvfdmJw5bW+qCpDgkZ3RUHUM/WE5a/JmOqkILwiJCZPM7W1v08ssaXyiyGCw+eFD5T1Sx20ELwgJSZiMOa2vqtygO9hgbU8f41r6oJGa0Wtk9EFyCApExzz+ffiCyXlcXbRV9F1zFn/lIS9W6HOxiZEgWT5CqD3q5jpTK3yGNG3qVpkNpQWREWnFotSgOruAHxF1fT/4b/nRo2ckkopQQ4oIQ4ich62aFEFIIseJ9fuctv+T9nU+7fRdIeyFIyTcIXWRyMr1ol10kUglEMDEOgEYwMY7HtJRSSCnPAngCKsEOAMwLIVag8mKUlnY73QWwsUGroQykdbQqs0iktSCehIpSDXSS4NzFEw6fGSmlLwhPSCnP9qwvFd/5TrqQZfRnKBe+o1VSyioSaQWiBo0kOF6A2qAYzEZ1S8rA7i7w5z8nq/vss7QaykyaZD1lFIm8Bikf88LgAwCklK951sMxTzy6EELMCyGWhRDLN27cyKmJepw5k3yEu90GfvWrbNtD8mdzE3jhhWR1yyYSaQWigZgkOB53xya8QU0/RP4tqNwYXbiaOGdqCvj442R1s3KeIm7w8stK8JNQpusgrUC8hc4NfjcJjhCi5hcQQsygM04BqLEKv7txFkApYtofPqxeH04CuxTVRIjk57YsIpFKIPokwflTT9HVnjrf96yIlTIkzrlzJ9mceJoLiJSHKotEak9KKeWbIcvOBb5fA/CDuDouox+8VAJQZ/3nPwd++Urnf1Jt/EBApmLh+kOErtYx6Kr8S1hEDQ1cxGW02wICErh4UUVzWVy02UTiCO22erFuZ8esnssiQVfrPuibgBI1NPA0liAXLnbEYWlJedi4evZJ5mxvJ6vnaneDAhGB2QkTuIjLwMKCEoWhIfV3YQG4fNnds0+sICXw3HPm9Vy8TJi8N4QkJ0pK7I9Im9YXm5SaF18EXnnFvF7et2S/5L20IHpIJQ4XL3avuHiR3YsB5uWXVXwQU1zK6EWBCJAkn4Rsy444LC0BTz2lLAe/u0GRGGheecX8nZs7d5Q7vwtwFsMjia9D+8VF4GJDjTPUakocAOAXv1DLALWc3YyBZmPDfAp0ZMSN5woFwsPUrPvWNyXE7YayEgAlCE8/DbzxhrIe/GUUB4JOlnIT92wXpj8pEDC/h8fGgP/5XwFIz0pYWuoIBWcuSAR7e+aXxeRkchf/LBj4MYif/cy8zt3w5kJ0uhI+FAfSB1OLYGOj2CC4Ay0QUgKvvmpep+sfzlwQQ0wvjyJnNQZaIEzD04eKg+8QxZkLYoCpVVCUUTqwYxA/+pFZ+X0nVAg1QxEcc+DMBdFkeBg4dQr4/HP9Opub+YcpHEhPyl6HxzjuvRf45JM+PxYUA0aGIQaYXio2bld6UvZg2rWIFAdg/xmmOBADTG/4H//YTjuiGDiBMJ0yctzAIhXAZDzi17/O95q0mjjHW78vSU5cHZuY5Ltstey1gxCf4WHgnnv0y6fN/WqC7cQ5QE+SHM06VnjmGbPyzHZF8uKzz8zK3/XFsYzVxDkevUlydOpkjpT7fZriyhOSJyYW68SEvXYEySNxTm+SHJ06mWNilq2uxpchJGtMLdY06QB1sd6biUuSE0bWiXNMzbGjR1NvkpBEmFiueVynVhPnRCTJ6VsHyD5xjok55sp7+GRw2djQL/v009aaAcB+4pywJDmhdWxhEmF4bAw4MLC+pcQVTLwll5bsjpdZTZwTliSnTx0rjI7ql81rZJiQOPb29Mt+8YW9duSROCdsfS6Jc0xGhW9FZRUlpACGhoCHHwbeey++7D332LMiKu1JaWKqTU/HlyEkT65c0S+bwVh+KJUVCBProciIPYREIQTwgx/ElwOAkyfttKGyAmESKSovpxNCTPntb/XLJkkwHUclBUJK4PXX9comTZVGSF7oumEnSdsQRyUFol7XLzsyYq8dhGSByYtcJj4UOlRSII5pOm9zWpOUBd1ByJdeyna7lRMIEwUdG7PXDkKyRHeWLet4EZUTCF0FNemGEFI0Q0PAuXPx5YBsfXoqJRBSKgXVoVaz2hRCMufdd/XKZfD60l0qJRC6VsGHH9ptByE2MMmPkdXsXKUE4rHH9Mp9+ct220EqTG8HP+fIQmtreuWyihVRGYFot4H3NV77GhriG5skIYuL3UmR/ORJi4u5NUE3pqrJ1Gg/KiMQurMXfCmLJEJK9VgOZk7zM6s1GrlaErpWhEmogygq8yzV7XNxcJIkIpg5reBs7rpatLOT3hGwMhaEzshtHjH8SIVxJJu7rkv1hQzCQVdCIHS7F0VmSSYVwKFs7jpd5StX1NhcGvJInDPvfS4Flu1LppOGrS29cnknPiUVwrFs7rrXctrXCawmzvHCyr3jRZCaCUS17kqmk5bjx+PL2AqoQQaEqGzuCwuFZHPXDZGfVrfSDlI+CeBt77ufBCc42Tjjfd701vvBap8IJNJJhe5IrUlsSkJCWVzszt7ui0RBCZubTWBysn+Z734XWF5Ovo20AlFDnyQ4PbEnZ6EiWgMqmQ4AzEopX0vTAF0TKu5AEqKFQ9ncdR6OV6+q3lDSfJ65DFJ6XY+3/QjWccl0TBLn6Oz4M88Ueh4JsYJu4hxdv4kwrCbOCXDBtxQikul0YZI4R6eP9eKL8WUIKSM6N3+a2TvbiXMghJgPiMMFhCfTSYyO49PwcJotEOIud+7El0mTLc5q4hxv+SUhxIoQoh6o05VMJ+n2b9/WK8fxB1JVdAbf0+TwtJo4xxtn2Ne8rBLn6Iwr6IoIIWVE1x9CZ8YjjFJ7UuqYVwU4uRGSKzp5XZK+uFVqgdBxI6VAkKpjMyN9qQVCx8X6yBH77SCkSHTe2NSxtsMotUDcd1//9fPz9H8g1UfHF8gkFWXXbyerVjw6Kc9ffdV+OwgpGp33Mh56KNlvl1YgdEj7qishZWFlpf/6eh3Y2zP/3dIKhM7gI1/QIoOCzhRmkuS+pRUIncCzSQdmCCkbOt7CScbjSisQOhF10niQEVImdGYykgSPKa1AxPGPfxTdAkLyQ2d8QTdkfpDSCkTcAUkyIENIWdF5aTHJa9+lFYiJiXTrCaka12ICOCZJx1dagYibxaCLNRk04t7J0Hlno5fSCkScjwN9IMigEZduL0k6vtIKBCGkm7g3NpO80Zk6HoQX+KWBiAC0Yevj6ugQ53+eNEgnIWUlbhoz92lOjbwY+9bH1dElzgmKTlJk0IjzhUiSpzPtc/ZJKEsA6OTFiFsfV0cLWhCEdBMXnDZJ8Nq0t1ENffJiRKyPq6OFjYNBSJmxMQbh5HNWJy8GHaUI6cbG1L/tvBhh6+PqaOXFiAuAkTRABiFlJS4uhG4+zyC282KErQ+tY4qNARlCykzc69y5v+4dlxcjbH2fOkbExXpgLAgyaBw+nG59GFbzYvRZnzovBl2tCelmYAYpdYh78STJiymElBkbvkGlFQhaEIR0Q4EIUOlByl51o9oRDWyMy1EgXGNxEbh4sSMKUqr/FxeLbBUpAS76QRRGXLoxm+nIrCEl0GgAS0sdkbh4Uf3faNCSIH158MH+6wt5m7MoKukoJQRw+bL6vrSkPgCwsKCWM01YsUjZfQ56/y8QnZufXYwAX/1qPu3InKBI+FAcisfxrt/GRnyZsTHz3y2tQJw6FV8myfvvheNfeEGCFybJnxJ0/XS61InG5aSUTn/OnTsnw7h9W0p1ZqI/m5uhVd2l3ZZyYUE1fmEh/P+82tHv/0EkeC78T57nJIb//Cf+frh1K7wugGUZcf+V1oLQ6XMlGZQpFCFU/PLgmMPly+r/Wi2fbobjpnRhON7103l7OUlTSztIqZN6zwHLz5zFxe7BL//CzONCDJrSgNqub0ovLDg1KJc7UV0/R0RC52GYKNNclGnhyieqi7G3F29SffSRkZVGpHTelC4EV7p+fYi7F1ZW+tWtYBdjaAj4xjf6l7n//nzaUjkcNqULwYWuXx90BuOTzGAAJe5iAMAf/xgf6391FZie7l+GeCwuAvX6/uWPPgq8+27hN0KhFNn1i0HnHYukIRhLa0EAeo4fDF6riZRKHN54Q32eekp9AODKFeDpp0s6qJMhvWLggDgAeqclqUCU2oLQUU4GjtFECOD115UYXLmiRALoiMTRo87cEKQbHQt5eDjZb+eROGfe+3pWSvmst+ySlPJZIcS8TBE8Rsc6GB/ng08bIVRXInhgX3+9s444x+6uXprJ8fFkv287cc4FAO94IjDj/Q8A80KIFai8GInR3elSvrhVBFFTecRZ1tbiy5w/n7yrbTtxzkxg2TV0gtU+IaU86wlLYnSj9DYasUWILw6+z0O7rf4G3YuJc+h0Hf7wh+S/n7aLUUOfJDg93YdZqIjWADArlMmaODenz+ZmvCXBJDoaRE3lAU5M5ZFwdJyf0sRGyWWQ0ut6vC29CNayk8T3MSHEhV5Lwhu3mAeAM2fO9P1tHRPr6FE+ALVweCqP7Ec37ur2tkU/iMAgY5Br/rgDYpLgeFyQ3Zm9IaX8b6/8TG9hz/J4EwDm5ub63toTE3F7oFhfB6am9MoONI5O5ZH9NJt65Wq15NuIFYiYWYa3AMx537sS50gpG973+YA4XIAai/AHJ88C+F2ilntMTuqVowVBqsZLL8WXCfN7M8Fq4hxv+SUhxIoQoh6o833PkliRCRPnBNE5CEeOpN0KIe4gJfCb38SXS2sAWk2c43VD9g2jpPF9CEM3nHerlbwvRohL3IrqzPeQdoC+Eo7Iun0sdqdJVXjoIb1ySR2kfCohEAcO6LlUc5CSVIG9PeCTT+LLffZZ+m1VQiAAvQN25w69Kkn5uX1br1xa6wGokEDoZi4uZTh8QgIcOxZfBkiWzbuXygjEyIjeGANnM0iZ0XWO0nEg1KEyAgHoH5T1dbvtIMQWuhawThBbHSolELpTOl/7mt12EGIL3cCzabwng1RKIHSDw1y/np3CEpIXuoOTui7YOlRKIAD9g/Poo3bbQUjW6FoFOgFkdKmcQOjGiHjvvWwPpFV6XyTRebEkSR3iLKur8WV8svT3qZxADA8Dp0/rldU12QolSaYrZseqHLpTm1lf05UTCAD4v//TK+d8OPxgpivdpLFJ6hCnMYmIloXvQxdRGXVc+URl1oojLtOQ//n61xP9fH4kyXTF7Fj5YjnZse61vLqa9PejM2sJ6fgTZW5uTi4vLxvXa7X0XU339hzPnyFldwPb7XivsCR1iDmLi+oR70fe8i22Wi2TLt32tv642s5Osrc3hRBXpZRzYetcvi1SYfJat+4JKISoSNP9hD1JHWJODt053WtzeNhS7NUo08KVT9IuhpRS1uv65tnubuLN2CNJ0tgSJJqtFBa7c82m/vW7uZl8O+jTxcgjcc6+JDlxdbLCxJvs4EEHH7BJIk0zOnW++Md3aamzLKNAv7rhFAGLgZCilEPnAxXK/nHv+zzUDd9bpg5gBSpwrVad4CeNBSGllDdu6Ktwq5VqU/ZIMghmeeCMeFiyIFZX9a/bZjPdLqCPBWE7cQ6wP0mOTp3MOH5cv6yz4eiSRJpmdGr7SGkt2ZDJFLxuZPckWE2c49GbJEenTqasrenPDzebZqYdGWAsdedMqm1tJdqENtYT58ieJDk6dUwS5+hg4no6NeXgWARxl4yTDW1smJW3nb3eauKciCQ5fet45bUT5+hSr+u/KutPZxOiRYbdORPrNY/wibYT54QlyVkOq2Mb0/fjNzbs9u0I6eXhh/XLTk6qYM22sZo4R4YkyelTxzom/TWOQ5A8abcBE4fhvKKiWU2c02d9polzdDHtr7GrQfJieFi/bJ4hEyvrah2F6Q3PXBrENqZDFnlatwMnEIDZSHGzqZ/ajxBTnn3WrPzOjp12RDGQAmGaUMTKSzBk4JESeM3gRQNrL2T1YSAFAjAPN0dHRJI1piEGirBkB1YghACee86szvPP22kLGTxMHzhFdXMrGzBGF9MTxbgrJC1DQ2aD5adPA59+aq89AxkwRhdTfXQ68hRxnhdeML/mbIpDHLzcAWxumpWnBUGSICXwy1+a1ykSCgSSveZNkSCmmFqfLmR/o0B4JFFqigTRxfRaOX/eje6sA01wB4oEsUGSa+Tdd7NvRxIoED1QJEiWJLk2ih53CEKBCCHJe/YUCdJL2cUBoECEcuBAslgQFAnik+RacDGZNAUigmYzWT2KBElyDbzwgpvXDgWiD0nNPRdPNMmHJOf+9Gng5Zezb0sWpBYIIcTjQogLQoifhKybFUJIIcSK9/mdt/yS9zcs3qVTUCSIDkld8A8fLtZTMo5UAiGEmAUAP4Ct/3+AaSmlkFKeBfAEgEve8nkhxAo68SqdhiJB+vGtb5lFhApy+3a2bckaq4lzAslyAGBGSukLQm8yHedJIxK2cxeQ4hAC+MtfktV1bcYijLQCUYNGEhwvQG1QDGajuiVe+XkhxLIQYvnGjRspm5gdSU/o2BitiSqS5pyWQRyA/AYpH/PC4ANQyXQ86+FYWDIdKeWbUso5KeXciRMncmqiHmlOLEWiGqR95b8s4gBYTpwT4O7YREQynVIRTKRkit/lsJ0Ridjh/HngypXk9cskDoD9xDkQQvjZtHzCkumUjjQicehQ5zdIeUhjNQwPlzP4sdXEOQFWe+p0JdNJ04YikTJd9i0h3PSeI93s7KQTh+efL6c4AAw5lwn33Qd88kny+ufPu/P2Hukm7bhRGUIUMuScZT7+ON0T4m9/UxdRWZ8yVaTVSn9jp+mGugIFIiOGh9OPKRw82BmfIMXgP/FNc6f04rhhrg0FImPSXhjb2+oC3d7Opj1En4cfTu4R6fPNb1ZHHIAMkveS/UgJnDmjuh5J8S0JTonap14Hpqfjy8Wxt+dGmLgsqdjuuMO//53NmMKhQ8qiyDsn4yCwuqqObRbiIGX1xAGgQFgli3EJn9FRdTEniXZFuvGnLY+FvhhgxuRktboUvVAgckBKs4zi/RgZoUWRlHpdHbusumy7u8D6eja/5SoUiJwYH8/2SeNbFKur8WUHmd1d4MaN7LoSQOdcHhiAEbwB2EW3kFI9dQ4fzub3fDO5XlfjFZwmVWxuAlNT2Xuq7u4OhjD40IIogKkpJRSPPprdbx492nmtfH3djaxMeSMlsLamjsHERLbisLU1OFZDEApEgfz1r3Zu5MOH1YV8+rS6Yao8sNlqqW7WD3+oZhGOHMl+G1IO7lTzgOmhe/ip4Le3s+8efPZZ9w1z/boadZ+aynY7edNoqBmirLppUdAHhQLhDKOjSih2duxdlKdPd76vr6un79SU++MWW1vKp+TOHdWVss3OjnJ7JxQI5xgZ6UyLTk7a206vFVGvK2tmdFT9LeoGabXUvo+Nqe95BRQTQg1sui6WeZOJQAghZqPiOnhxHxoAZqWUr0UtI91MTCihaDaBb38beN9y1IywJ/OHH6obdXRUmfRbW2qKL00MDED9TrOpftvvXq2tATMFxBb74gugVqPFEEVqgfCCxVwCcC5k3d2w+EKImWBY/OCyMgeNsc3kJHD1qhqRX19XF3NePPhg9LqPPlKDn0NDnUHQ3V015tFqqZt+ako9mYVQ4d1rNRU7wwVu3szGk7LqpBYI70aPctd5EsDb3nc/LP6xkGUUiBj8EXop1Q2Y9nXktNx/f7HbT8LkpBKGQR94NMH2NGcN+8Pihy0jBvimeaulvARJNBMT6hj5btEUBzM4SFlifM/JoIOQ7am/snDzprK4Bs2xKWvShr2Po4HwsPh9Q+V725wHgDNnzmhshvjCIKXyExgaAk6dGpysXo88Arz9ttrvtIOopEPasPehBMLeh4bFj1jWu803ARW01nT7g44/kNlqKd+B3V3lsVl2B6kg588Dv/+9EsQjR6oZi8EFspjFeBzAnBDicS8ZDqDC3p+TUr4vhJjrDYsftozY4cCBjpntT5vu7XUCnNhwTbZBs6kEb3RUWQgUhHxg2HuCnR01gOePZ2xuqk+/aU4bfPqp8kcYH1fTuu22chwbG8u3HYNGv7D3HMIhGBnp9gnwPTh9J6aNDfXEHhlR4xsbG+pJnsQPYm9Pldna6tz44+N0VHIVCgTpy+ho99RgUv8LVxykiBnsyRFCIqFAEEIioUAQQiKhQBBCInF+mlMIcQPARxpFjwO4abk5tuE+uEMV9kN3H+6XUoZG3nBeIHQRQixHzeWWBe6DO1RhP7LYB3YxCCGRUCAIIZFUSSCMXypzEO6DO1RhP1LvQ2XGIAgh2VMJCyIY6zJk3eNCiAtCiJ/k2SZC8sbGfVB6gfBeG/+viHV3g+YCaPQ7gEUTdwKFEJe8v2EBfApBo83Oi3MZj3sYtu6D0guEt9P9guY2vO9+gFzn0DyB80KIFaj9KJy4NpdBnMt43KOwdR+UXiBiqKEcAXJ1TuATUsqzmqH+8iCuzWUQ5zIe9yTUkPA+4OveblBD/AmcFUIA7iQbqqF/m+PWu0AN5TvuueK8QFgKmls6AlnJHhNCXCj5E600VOS4N5DwPnBeICwFzc2dGKFroM8J9OJ+wov5eQtqX4qmgf4XXdx6F2igfMddmyzug9KPQQSD5gYW/wkAAkFyCw+QK6V8M+Tjn6i30Ln47p5AIUTNW3YNnZN6FoALQTrj2hy63jHKeNxDsXUf0FHKETwL4xqAGd9q8oKJngusX/XWO9EX1mxz13rXKONxzxMKBCEkktJ3MQgh9qBAEEIioUAQQiKhQBBCIqFAEEIioUAQQiKhQBBCIvl/CU/DjFCw9bgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "uepd.plot_object(X, X_noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgm_tot = uepd.alphacomplex(X_tot)\n",
    "dgm_true = uepd.alphacomplex(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_points_samples = 100\n",
    "nb_repeat = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resolution parameters\n",
    "m=-0.1\n",
    "M=1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "diags = uepd.expected_dgm(X_tot, nb_points_samples, nb_repeat)\n",
    "h = uepd.tohist(diags, m=m, M=M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAF5CAYAAABdm6YkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABASElEQVR4nO3dXWwc6X3n+98jcWzNi60WNYqd2HGkVoLEiwziaVHnHBgHOGulJ77JTTaUtYAv1rB3qEtBgVeK4Rte2JkRsxhBN96Ic5K7uVCGyf2xOGMsFjmbA1E9OjiLeNdrNqUZzZs04rQkSkNSJJ9z8TzFri5Wv7K7uqr6+wEaXV1V3V2U2E/9+PS/nsdYawUAAACge3uGfQAAAABAVhGmAQAAgB4RpgEAAIAejQ37AAAAANLEGFOSNC5J1tr5IR8OUo4wDQBADhhjLkgqSyr5VfOSaqFdCpKK/iZJJ621c0kdX8aclvQdSVVJxzp5Qpt//4K/r0k6b62tdvnccUnLki4Q7tPHMJoHBq0fjYQxpiDput9ndoCHCwCZZoz5VNKCtfalJtuLkq5KumytnUn04PrMGFOy1lYG9NrnJJ2y1nYUpkPPa/rvb4wpS3pT0itx//Ztnlv0z61aa092c0wYLHqmMXDW2vNSZ42EMaZVIzGu+l/3AIB4y602WmurxpjTkmLDdsackjSQMK3GXv1uLDd7rrV23hhzXtJlY0wlpgOp1XOrko4ZY64bY94kUKcHFyAiSS0bCf/Xf9EY82bM9pq19kDWe1EAIA3yUCrgO2EKwz6OboW+XT3d40v8saRJY8xknw4Ju0SYRtrQSABAMu4N+wB2qdcwmhal9rvsZK2tSZqTdKGvR4OeEaaRKjQSADAYxpipyKpq7I4pZ4wp+J/l3LCPpRf+GiDJXT/Uqyty3+SWd39E2C1qppFGV+R6p8vBV5H+67yipOVBXWwCADl3NPwgGMnDGHNZ7iLxYJSPY0E7678lDErvKnLfHo77dUVJs5IuSwq+TTzoX/t83AH4i/pqwftIejOu5MSH5aOq955XQyOPTEk67l+nHC4NjNYRd/F+k/7nqcmVjgzyPBP8UbObTqPg+F7S7kI5+oAwjTSKaySCxn5OUrSxLMg1TjW5RvCg3JXStegL+31/JGnR71uTtCBpQvXauznFnyiOSqqFLqhs+76hq6+Lkv5e0nm54ZYkP9yStfa03zdoYAtyJ4qXY16rrPrFReP+fpnhrQC0Emo/dgi1QVcljYc7LKy1c8aYqqSXQsO51eQvhJMrVSiFr2cxxpwzxizKhfJaaP11uXatElq3aIw5HQ64PhxXw4HcGFM0xkxaa+eC9/L7LQfHH/Mzd/p+lxVq2/26klx73Ve+nT+txn/PXgTngWLLvZAIwjTSaEcjYa19Ke7CRG8q0pCXJV03xhyLCdTXJW03pL4Rfclae9KH43L4imn/nLK1dsYP8TepegPb9n1Dr3W1yXPe9K+7KOnvg+f53pS35AO3P7bz0ZOG/ze50uTfBcDomvDt27hcW1pS+7KOk5KWjDFTwUVyvi1qFvyqkorRP+Z9e3laruc1COoX5EZzivb4XpDrsDjq95uUa3MPRPY7r3qHSls9vF+0177i/4gYV29K/hgCBX9/PfpevbDW1owx4dfFEFEzjdQJBeBCZNOOxtwH2FOhGrTwVerfiew7JdfrEv5K7LJcSUnBjxgSbqircuH57/3rnlc93Hb8vqHX+o5cT3fYVbne7Wok+M/JNcbB65cVPxLKKzHrAGDBWnvaWnvSj5R0VG2GevNt0Em5YduKvne22qYHtdm2y6qXM0iuvvlqzH4LcrW/Bf/4gnybG3FdHQbpHt6v2esudvF+UVVr7fnQ7bS/9WWehPC5B8NHzzRSp8tGYlmu12VcjSeKiiL1gXJBONrwB48nFF93VguH3NByN++7/V4xPeXLkgpNxhpV6PWrkl43xlyJfAVb8b0TANCUH1u6bW2tHwd5Vp1P6tJsTOug5jo8YsXxmPZqXI1htigXnKPH1XEIDb1np++3m9A8LEGPeSYvIs0bwjTSqONGwgfLA9J2b3FJLnwWtbORr2lnb/d23XGTt4g9hi7ft+VrqYOJAXxonpcrI6nKnQyuWmvnuSATQId29NSGSzoC1trTfpKtg7t4r6C9C9f0XmnSXgUlJcG+LSed6UDwOkm93zAExx7X+46EUeaBNOqqkTDGTPn65oK1dsafGOKC61XtvFijKNf73CyQNm1ku3jfvvBXqQcXZU5KuupnwioM6j0B5Ef0G7CgjCO6ny+Je1nSuV0MvRa0tVXVLypvWX8cKifp6aI6P2ReKan3G7KTcucuLj5PAcI00qjjRsI3+pcl/XG7/f2JpOIvqAmclztpdKWb9+2H4GtL3xN92l/AEpSTMCY3gF6UFQnTvrc2GCHotKQ32/zB3iywBm1WJRRaY6cvj7x+Vc1L5TpRTPj9EhcaSYprZlKCMI1U6aGROC1pNqYWuRB6zanQa78iqep7lafkRsjoJQx3/L59MhHtIfInjJNqMtwVALRxOubiwsmgTfTfti1Ier3FazTr1T0tKVxvfV71saijwm3lBcVfxK1IR4hUv3YlMK76HwfdvF+zNvS40jlaxpuSKh3UsyMhhGmkTbeNREnxF49MaGePyYR8j7e1dtbfeq037uZ9++VkdIU/EdYG9H4AsqltG+SHzYuuu6qdF2KflhvxqGnngB9eLvz4giLjNvs2vRp930hZRhDg55vsFw3+b8q1uYFy0KZ3+X7RbyyDa2GCi8y7Nd7j84LnFuI2+BFWrkuSH6EFKcEFiEhSy0ZCrmFcbtFIFGOePy/3VV547OagZzv6FV9V7ivLk/K10HETu4QU1Lxmupv3lZr/7K3WK7LtO8aYC+GeJN/g7zgpAhg9PsSW5dqNhpkBvYJcOxr05gbj7b8Zet4pNc7+F/TaXvZt5+XIt3kLcqF1Uq7dOirpXlw77ucLOOcD7qL8SEXRWm4/7v9UaD/JdbJEx7OeN8a84o//WuS4u32/c6ZxtsSq3Bj+F/zFmC+3+xYzGJNb7t+x5P84ebOTkUg6/L+T3MRg1EmnjLHWDvsYkHOhRiIYrijaEBT8fbSRDp5fkPuaMej9mJd0MjTBSRAmg+GUFvzoF1flGsTLtj41bnTsU8k1nvNyJR9VH+yDGReD470SPbZO3te/djDZS3Dsp/37hE9gFflG0jfIp+VOeMHIHdf8crinpCDX+9OXcUsBoBumzQyEwKggTGNk+FB/LRyKfVAvyoXX70g60qa3GgAgwjQQoGYaI8F/BVmI+Zqw5q82Py3Xaxx74QsAAEAcaqYxKqIzFcZZVjYH7weAYSgonaNdAImiZxojwdcVl5pNQBBcPMOFHQDQmh9V4qrcNR9lY8zV0GyCwMihZhojxY+4cVSNw9odVaSWGgAAoBOEaQAAAKBHlHkAAAAAPcr0BYjPP/+8PXz48LAPAwC6srW1pXfeeecTa+2hYR9LkmizAWTVzZs39cknn5i4bZkO04cPH9bCwsKwDwMAOraysqJbt27pD//wD28N+1iSRpsNIIsePnyob33rW023U+YBAAkJgjTXqgBANjx8+FDvvvtuy30I0wCQAII0AGRLEKTbtduEaQAYMII0AGRLp0FaynjNNIBkPXnyRLdv39bq6uqwDyUznnrqKf69ACBDugnSEmEaQBdu376tL3zhCzp8+LCMib2oGSEbGxu6c+eOHj9+rD17+CIQANLuwYMHeu+997r6JpHWHUDHVldXdfDgQYJ0BzY3N/XkyRMVCgXKOwAgA3oJ0hJhGkCX8hykz58/r5mZmV2/zubmptbX1yXl+98LAPKi1yAtEaYBZMzc3NzAXvvUqVO7fo1wkAYApN9ugrREmAYwQG+8IR0+LO3Z4+7feGN3r1er1XT16tV+HFqsQqGwq+cTpAEgW3YbpCUuQAQwIG+8IU1NSY8fu8e3brnHkvTd7/b2mgsLC1pYWNDc3JwmJyc1NzenK1eu6NSpU6pWqyqVSrpw4YKuXr2q8+fPS5IuXLggSZqZmVGpVFK1WtVUcCBesK1SqWyvq9Vqmp2d3V5fKBRULBZ14cIFnT9/XpVKRcViUYVCQW+++ab+6q/+Ss8880xvPxgAIHH9CNISPdMABuTHP64H6cDjx259r8rlsorFoiYnJyVJk5OTqlQqmpyc1Llz51Qul7f3PX369Pby+fPnVSqVVC6Xtbi42PCaQWAul8sNz19YWFChUFC5XNa1a9c0NTWlcrms5eVllctlTU5O6sqVKyqXy3rxxRf1z//8z73/YACARPUrSEuEaQAD0mz21TazsnatVCo1PI4r1ahWq6rVaqpUKjp69GjDtuvXr6tYLO54ThCs5+fnt3u3JWl8fHx7uVgsanNzU1/84hd38yMAABLUzyAtUeYBYEC+9jVX2hG3vh8qlcqOIB1WrVa3l48fP65isahSqbQjOB87dkzVanXH+kqlst0THhYO61tbW9RIj5SNNsvNtsc9lhpPwWNN7pvtH4dTOtBOv4O0xCcPwID89KeNNdOS9Mwzbv1uFItFzc3NqVwua35+XpVKpSFYHz9+XPPz86rVapqfn1e1WtW5c+c0MzOj5eVlSWoo55iamtoeDq9Sqejq1auamppSqVTSsWPHVCwWVSwW9aMf/UgLCwvb7/fzn/9clUpFS0tLmpub04EDB3TixInd/XBIqQ1Jq6H71ZjHa5H9os8Jr5Ok5yTti7l9PrQ8Fropct8sfLe7B0bX/fv3dfv27b6P/W+yPJnAxMSEXVhYGPZhACPjl7/8pb7+9a93vP8bb7ga6XffdT3SP/1p7xcfJm1mZkZTU1MqFAqqVqu6fPnydrlHt6N2/OpXv9LevXsb1r3wwgvXrbUTfT3olMtum70qacXfav7+UWhdsH7NL69G7h9F1klSQdJ+uVBdCN2eDS0/p3pw3qvGcN3LDRhduw3S3/ve97SwsBA7cQCfLgAD893vZic8R5VKJc3Pz6tQKKhWq22PQc3wd6NoQ/XAfM/ff9Lk8X01hO6Nx9JjuQwdZGpJ+qJ8lt4j6Xl/+1JkuaB6D3Xc/V419mC3WpY45WNUDapHOsAnCwBihEtBAgTpURWUZ9TkQvNHkj729+HlT6T1demBpIdyIfqhvz0IPd4j6aCkcUmFLWn8jnTwjvTFX0v6qr/V5EJ1XPlHs/KQDe0Mz/tC64HRM+ggLfHpAoCOEKRHWdAzfV8uTL8v6Xbo/qb06Va9Y7rmbw8Vv26vpEP+9iV//2VJX1qXDlWl52qql4/sU72+Orj/fMy6fXJlJs+GjjsI0sGNUz5GSxJBWhrSJ8sYU7LWVppsm5RrckrW2plEDwwAYhCkR124Zzrohb4l6aa0/kG9g3rZ3+75+1rkcbC8V9Jv+dtXJP12aP+Hkn5rWXo+uLDxObmA/FzktiYXqldD6+JGFRlT/EgiQL4lFaSlIYRpY0xZ0gVJx2K2lSTJWjtvjCm2Ct0AkASCNOoFzzW5nmnfG/3ZB9IHqndQ321+21yuP9wj6fDH0rPvyQXpu3Ih+6F/myeS9Fh6/teqX4gYXKj4bGg5CNDB7bnQMQf10quq91ADoyHJIC0NYdIWa+283N/gcU7JtVaSVJW0s2gxQW+8IR0+LO3Z4+7feGOYRwOgX+bn5/XSSy+13Of8+fN69dVXCdJQvWd6RdtlHhu+R/o9uU7qqqT/LulfJP1/km5IqkgP/of0q2Xp/5b0XyT9Z3/7L5J+9Ylk3/H7/jf//Kp/zQ8kfbourd+RdFMurd+W6xn/JHS7p3odSXh0kfDwfQRpjI6kg7SUvhkQC2oM2geHdBx64w03Ru6tW5K17n5qikANdCXamPWhcZubm9v1PuVyOXamxLDJyUltbm52c2jIrSBMBzXTH7ne5A/87Zak/ynpf2g7TD9Ykv77Y2lB0n9VPUi/7e//q6T/R277gyX/vF/Lhelb/nWD8pFPt6T1ZTUG6iBU1/wt6DkPgvSaCNMYNcMI0lL6wnRq/PjHjZNNSO7xj388nOMBMmd6Wjp7th6grXWPp6d7fslaraarV6/uep92Njc39eyzz7bfESMi3DN9X9p4UK+HDgLvB3I9yjelOysuE/9PuYz8L3Idz/9N0v8rqSLXGf0vqndGr3/gnqsP5IL6shovXHwg6bN1NQ7DFw7Q4d7oaK80gRr5N6wgLaXv0t6a3GBBkuulvhfdwRgzJWlKkr7Wr3mJY7z7bnfrAYRYK9Vq0qVL7vHFiy5IX7oknTnjtpvYse9bWlhY0MLCgubm5jQ5OSlJmp2d1cTEhBYWFjQ1NRW7T9BTvby8rKmpqZbv8eqrr+qFF17QjRs3ttfVajX93d/9nb7xjW/oxo0b2r9/v44cOaLXXntNf/EXf6EbN27oyJEj2r9/v/7xH/9RP/nJT9r2fCNLgtEyCpIOSmPj0qFlV+O8Klfj/ETSU5K+IP3GB9LBj93J7JDcdYa/o3pH9h5JfyDpd/3970na+7uSft8/+Kp/0iFJz7jX1LOSzBfVOMHLQTVO8BLcwkPmMWkL8q9Wq+n9998fSpCWUtIzbYwp+MUrkop+uShpPrqvtXbWWjthrZ04dOjQwI6pWU4fYH4H8sMYF6DPnHEBes+eepC+eLGnIC258oxisbgdkmdmZlQsFlUqlVQsFjU7O7tjn0qlomq1qsnJSV2+fLnl6//N3/yNXnjhBZ04cULf+ta3ttdXKhXt379fJ06c0PXr1/WDH/xAJ06c0PLysk6cOKE/+7M/09zcnE6cOKEXX3xRlQrXTedLMF7zs3JjP3/Z5dgvq56Ui3JB+A8k/Stp7x9KX3teOi7pf5X0TUn/u6T/Q9IJv/xNSX/wjNtX/8o/tyh3UeKX/O3Lkp57RjK/JZeyvxza+LzqMykGsykGt8+rcfIWIJ+GHaSlIYRpP/TdhL8PvCVJwcgdfsSP2jBH8vjpT6Vnnmlc98wzbj2ADgSBOmwXQTrOtWvXVCy6v7+LxWJseUepVNqezXB8fHzH9sDm5qauXbumI0eO7Nh24sQJSdLbb7+tn/zkJ9vrw68XPK/VeyCrxlQfQeN5SV+RzNfqYfqwXAj+fblQ/IKkP5JUkvZOSF/5Hekbz7kgXZYL0//bU9JvHJFL2iX/nN/3r3NYLlA/N+7eR4flgvRX/O3Lqs+UGO6dLoieaYySNARpaQifLmvtnKS5yLpjoeXZpI8pTjAF8o9/7Eo7vvY1F6SzOjUykLigRjrs7Nm+BepKpaLjx4+rWq2qWCyqWq3q+PHjO/ZZWFiQJE1NTenChQuq1Wo7SjCC4e9efPFFLS0t7QjU77zzjk6cOLFj/f79+3f9cyALxuR6eguqd0mvSZ/bkH77A+lzciUeBblixS+pcZzpmrs9vSx9pab6ONNBGA+Wvxpa99xv+BXhcaaDSVri1oXvn9XOMA3kS1qCtMQnrKXvfpfwDPQkCNLh0o7gsbSrQF0sFjU3N6dyuaxz585pZsbN7VSpVHTu3Lkd+ywvL6tSqahSqWz3UBeLxe11f/RHf7Q9/N0PfvADvfbaa5KkGzdu6O2339b3v/99vfjii/rmN7+pI0eO6MiRI/rhD3+oSqWiGzdu6J133tEvfvEL3bhxQ0tLS5qbm9OBAwe2e7ORB0HP9HNyvcFr2r6gb2yf9JXb0jPrrkj6YeRWi3kczIAYzH4YvY0F3d1f0c5ZDj8feRzdFtcrDeRLmoK0JJm0HEgvJiYmbNDrBGDwfvnLX+rrX/96ZztPT7uLEIPgHATsQmFXI3r0U6cTsrz22mv6/ve/r0KhoKWlJf3t3/5tQ7lHO7/61a+0d+/ehnUvvPDCdWvtRNcHnWHZbbNrqg/b8XHkPrz8ibS+7gbZeKzm93vkerHHI/djvyHXNR0u5YgLyeGwHN4+FrMtvEywRvYNK0h/73vf08LCQmwvEJ8sAIMxPd04akdQQ93Hmund6GZmw2984xv6xS9+of379+v+/fv68z//8wEfHdIlCKMFScHY40E5RUH1AHxf+tyq9LkV6UBoApWN9fpodWtyYfoLkp7bo/pFhEH9c7QWOhyEw73N0QsMo8vRx0D2pa1HOsAnDMDgRINzBoO0JEo2Rl5Q5hGM1RzUTxfkgvSKXO91eKKU0HjPYxvSc6vupjX/ekEQD4/AER2R4zm5mpCxJjd1uQ3IrrQGaYlPF4AR022QBhp7poM65VW53uPobINS4yQpceuk+FKNuNKM4P0VWW5332wdkD1pDtISnzAAXbLWyqSkh7lbwwjSaW380Y1wyG0XlOPCbavAC6CVTz/9VB988EGq29JUTNoCIBv27dune/fupbpRa2ZYQbpWq2X2jw+EtesdjrtFa5gptwC6kYUgLfGpBtCFr371q7p9+7bu3r077EPpytbWljY3N4fSIBtjCNO5MibXCx13H94HwG5kJUhLfOIBdOGpp56KnSEwzVZWVnTr1i3t2cMXcei3TgI1p1mgW1kK0hJlHgByLAjSWWmQkXadXuTXah2AVrIWpCU+3QByiiCNwQp6o4PluO0AupHFIC3RMw0ghwjSGJxob3MnQZpgDbST1SAt8QkHkDMEaQxeuFc6eNxqXwCtZDlIS3zKAeQIQRrJCZ8+N1psA9DK8vKyPvzww0y323ziAeQCQRrDw6kU6EUegrREzTSAHCBIA0C25CVIS4RpABlHkAaAbMlTkJb4bgpAhhGkkR7Ruul+4TSNfMlbkJb4lALIKII0hm9QAbrZe3DKRrblMUhLfDIBZBBBGsMVF6KT6JluNUkMkG55DdISNdMAMoYgjeGKhuaNmHX9fr+49wSyI89BWuLPWwAZQpBGumw0We6n4DS9IU7ZyKJ79+7po48+ynW7Tc80gEwgSGP44sJzkj3TSYR3oH9GIUhL/JkLIAMI0kiXJMNttGeaHmpkw6gEaYlPJICUI0gj/eLCdT8QnpFNoxSkJco8AKQYQRrpEQ3Mze77/Z5xZR6DfE9gd0YtSEuEaQApRZBGugyrd5heaWTHKAZpiTANIIUI0kivsTb3wXKr/Vo9N7xPu/ciaCM9RjVIS3wSAaQMQRrZEa5pjqttbhe4m21r9RpA+nzyySf6+OOPR7bd5tMJIDUI0siGIDyHlwd5Oo0L4UA6jHqQlijzAJASBGmkX7NQO4wgTajG8BGkHT6NAIaOII3siOuVDpYH/b5AehCk6/h0AhgqgjSyp1mIHsRQdc3qsIHhIUg34lMJYGgI0siu8MyE0XWDfk9geO7evas7d+7QbofwyQQwFARp5AOnUYwOgnQ8LkAEkDiCNABkC0G6OcI0gEQRpAEgWwjSrRGmASSGIA0A2UKQbo9iLwCJIEgjn5qN4MHpFdlHkO4Mn3YAA0eQRv60GwYvqfGngcEgSHeOTzmAgSJII3+iQTr6ODr2NKdaZMudO3d09+5d2u0OUTMNYGAI0sifjcjyRpv10W1AuhGku5f4n8vGmElJNUkla+1Mi+1Fa+1sskcHoF8I0sifZgE5GqjDE7rQK43sIEj3JtGeaWNMSZKstfOSasHjyPaq316NbgeQDQRp5FuzIB3XSx33HCB9CNK9S7rM45Rcr7MkVSWVY/a54O+L1tpKEgcFoH8I0hgdzUo8mm0D0okgvTtJh+mCpOXQ44PhjT48V40xi5H9AGQAQRr51ywsx/VME6SRfgTp3UvVBYjGmIJcz/VlSa8bY4pDPSAAHSNII992E4wJ10gngnR/JB2ma5LG/XJB0r3I9ilJr/gLE09Kmoy+gDFmyhizYIxZuHv37gAPFUCnCNLIv7E299FlIN0I0v2TdJi+IinobS5Kmpe2e6QbBBcpxqyftdZOWGsnDh06NLgjBdARgjQwpvahOloOEh0+D0jOxx9/TJDuo0T/jLbWVowxE8aYsqRa6ALDtyQds9bOGGPOGWOqksYZGg9IN4I0RtOY6sPedTPTYdy+zJSIZH388cf65JNPaLf7KPFPb1xAttYeCy3vGHsaQPoQpDF6ouE5PI50N73MzUI1gRqDRZAeDD65ALpGkAakznumm003ziQvSA5BenBSNZoHgPQjSGO0RcNutF46qtkELkw3juQQpAeLMA2gYwRpQIoP0O1CdYCZEpEsgvTgEaYBdIQgDUSNqbORPALMlIhkEaSTQZgG0BZBGmgnHKo7GRovug7or48++oggnRDCNICWCNJANzq9gJASDwzORx99pHv37tFuJ4QwDaApgjTQq1YzJUYft+rRBrpDkE4en1wAsQjSwCAw/TgGw1qrjz/+mCA9BHySAexAkAb6ITpTYtwpl15p7B5Bergo8wDQgCAN7FazEo5mpR70VqN31lpKO4aMTy2AbQRpYBBazZTIaRi9C4L08vIy7fYQ0TMNQBJBGuivTkIzwRq9I0inB59cAARpYCCCU2yzXum4fYH2CNLpwqcXGHEEaWDQoqE6uh7oHEE6ffgkAyOMIA0kiVMudocgnU7UTAMjiiANANlBkE4vwjQwggjSAJAdBOl0I0wDI4YgDQDZQZBOPwq4gBFCkAaA7LDW6sMPP9Snn35Ku51i9EwDI4IgDQDZQZDODsI0MAII0gCQHQTpbCFMAzlHkAaA7CBIZw9hGsgxgjQAZAdBOpsI00BOEaQBIDsI0tlFmAZyiCANANlBkM42hsYDcoYgDQDZYa3VBx98oFqtRrudUYRpIEcI0kAabbTYxmk4H1r9H4c1/n8TpPOBTzGQEwRpYJiahano+g01nnqjjzktZ8dGzHKz34Od/6/W7iVI5wSfWiAHCNLAMMQFp7iA1Ww5egoe085wjXSKhueNJstjavx/dfcuSN9WrfaAdjsH+MQCGUeQBpLWaYiO660Mh6y44Eygzo4NxYfouDBd//90pR13VKs9pN3OCT6tQIYRpIGkxZVtRJfb9VpK0YBVD9ESp+YsiQboVmHa9Ui///6nun//Me12jvCJBTKKIA0MW7Mg3W65GU7J2dGsR3pV8WF6n6y1ev/9h7p/f512O2f45AIZRJAGhqHdBWetvvaXmoWseJye0y/aE70aul/1++yT65G2ev99o/v3N0SznT98WoGMIUgDw9ZtkI77+j8I0auhZU7J2RD3/x4N0mvbe1i7T++/v0f37+8hSOcUn1wgQwjSQBrF1VG3C9PhEB1XT81FiOkW9/8aDtKuZ9pa6f33n9b9+4YgnWN8UoGMIEgDw9Tq4sLgvlnAatYzHQTq8GtwWs6W6P9tEKRXfZA+qPv3n5K1ZpgHiQHjUwtkAEEaGLbwaBvdPi8QHa2j2WQthOr0aVXSox2PXZD+su7ff0bW7knmEDE0fFqBlCNIA2mycwKOzsNvuBd6n+LGIEYaxF1oGl6O+8ah/v9Y75EmSI8KPsFAihGkgTTrNlBHx5IOB+loqI7rBeeUPXidjBEeLeEJjMla6fbtZ/TgwRhBeoTwyQRSiiANpFU0FG/ErIvbfyOyrl3PdDSgN3t99Fe7kVmknf+fkrVjun17TA8e7KVGesTwiQRSiCANpFFciA4vx51SWz0nuI8L1K0COjXVg9duNJZokJZu396jBw8MQXoE8WkEUoYgDaRZu3AsxYfduMAdvRAxuI8rHSFAD16rHum4GmnHBekxPXgghr8bUYl/Mo0xk5JqkkrW2pmY7SVJRUmy1s4le3TAcBGkgSzoNEQ3q6luFrbD4nqmCdfJaDa8YThQu28TXGnHHoL0iEu0Ot4HZVlr5yXVgscRp32ILjbZDuQSQRrIkrjSjGgddPi2r8lyu5rp6HIvw/Ohc80uMgwmZHkkN5b0hqzd0O3b8qUdQzlYpETSf9aeknTVL1cllSVVgo2+13pRkuJ6rYG8IkgDWdXuAsFWvdhRzdbRA52suJ7pnTMbulE7uNgQCfdMSypIWg49PhjZflzSQWNMyRhzLrGjaiYabAg6GACCNJAXrXqlw7do73T4+YFmQ7LF3aO/4sJ0eGbDVd2+/awePGBmQzhpHATxnrW2Im33VDcwxkwZYxaMMQt3794d3FFMT0tnz9YDtLXu8fT04N4TI4cgDYyiuMCNdGj9B4rrkf5NPXjwDEEa25IO0zVJ4365IOleZPuiXPmH/P3x6AtYa2ettRPW2olDhw4N5iitlWo16dKleqA+e9Y9rtXooUZfEKQB1MWN7NFstI+46cjRH82+Wfi8rN2n9977qh48eI4JWdCgq0+iMeZluZE2rlhrbxhj/pN/XJF02Vp7s81LXJE04ZeLkub96xastTX/eDK0/Vo3x9c3xkgXL7rlS5fcTZLOnHHrDX+NYncI0siiPpwDECtu1A8maBmenYHaBeln9PAhMxtip25/I5Ylveob0ZclFa2137bW/khS25E3QuUbZUm14LGkt/z2qtwoH5P+8fCGxgsH6gBBGn1AkEaG7eocgKhmk7wE963GpMZgRP/998naz+u99z6nhw+fIkgjVrefypq19r5fnpR0ObTtfsz+O1hrZ2PWHWu1fSiC0o6ws2cJ1NgVgjQybtfnADTTyWQwzZbRX/VvCKy1eu896eHDPVR4oqlu/8QK1zi/JF+m4eXn1yxcI33mjLS15e7DNdRAlwjSyIHROAckqlVYJjwnq7FO3dq9eu+9z/Tw4SanfbTU7afzqHG9sv9W0py19oEkGWP+jdzFhflgjFQoNNZIByUfhQI90+gaQRo5MRrngMTF9US3qpkmWA9W0CN9Tw8fPiFIo62uPpHW2n/wdXKL1tq/NMbslzQlN170ryW9PYBjHI7padcDHQTnIFATpNElgjTyYqTOAYkLB+rgcdw+GDRr9+rddz/Uysoq7TY60vUn01r7emj5vqS/lrZ7JvIlGpwJ0ugSQRp5M1LngMQRlofNBelbWll5TLuNjnX9yTXGnJAbCimsIDdV+D/24ZiAXCBII484ByCvrLV69913tbLyiHYbXel2nOlX5RrRaszmQj8OCMgDgjTyiHMA8qoepFdot9G1bnumr1lr/zJugzFmOBOsAClDkEaOcQ5A7hCksVvdDo1Xa7bBWvsPuzsUIPsI0si5WrMNnAOQRQRp9EPX40wbYw7HbTDG/HD3hwNkF0EaI4BzAHKDII1+aVrmYYz5D5LK0dWSisaYA2qsmTOSXpT0H/t+hEAGEKSRN5wDkGcEafRTq5rpg5Lm1DjDVTMHJMXW0QF5R5BGTnEOQC5Za3Xr1i09esSoHeiPVmH6irX2nQ5fZ8kY80o/DgjIEoI0coxzAHKHII1BaBqm4xpRY8wXg+lj/eMjkkqSPrXWMvMVRgpBGnnGOQB5Q5DGoHR7AeJU+IG1dsla+w/W2reZ/QqjhCCNEcU5AJlEkMYgdRumW4nOiAXkEkEaiMU5AKlEkMagtZy0xX+FV5b0kqT9cldxvxTdTdIRSZcHcoRAihCkMUo4ByDrCNJIQsswba1dkvS6pNeNMX8vd1X3bMyuVWvt/QEcH5AaBGmMGs4ByDKCNJLSzXTiL0sqRy9K8T0XByTRkCK3CNIA5wBkx9bWlt59912CNBLRcc20tfZ+k+li70kyxph/37/DAtKDIA1wDkB2bG1t0SONRHXTMy1J8lPJliUVIptekvR/7v6QgPQgSAONOAcgzYIg/fjxY9ptJKar0TyMMX8uVy/3u5L+F7kLT573yyf7fnTAEBGkgUacA5BmBGkMS7c900Vr7Z9IkjHmj621bwUbjDEnJDFoP3KBIA3E4hyAVCJIY5i6HWe6Glou9PE4gNQgSANNcQ5A6hCkMWzdhulxSfIXmswbY66EtkXHHgUyhyANtMQ5AKlCkEYadBWmrbWvG2NelvQdSVbSm8aYLWPMPbkruoHMIkgDrXEOQJoQpJEWXY/mYa19XW4Qf0maM8YckDTuB/cHMokgDXSGcwDSYGtrSzdv3tRnn31Gu42h67bMQ8aYl40x14wxr4RWH+njMQGJIkgDneMcgGEjSCNtuh0a7z/IXXTyHblpZYOB/N82xvyb/h8eMFgEaaBznAMwbARppFG3ZR7VYAYsP4VsmOnPIQHJIEgDXeMcgKEhSCOtui3z+DS0HG04D+zyWDCKog1iQg0kQRroCecADAVBGmnWbZg+5qeSldyV3JK2B+unIUV3pqels2frAdpa93h6eqBvS5AGesY5AIkjSCPtuirzsNb+tTHm5/7rvZoxpiqpJKlirT01kCNEPlkr1WrSpUvu8cWLLkhfuiSdOeO2m/5/a0yQBnrHOQBJI0gjC3oZGu9PjDElScf8qlette/097CQe8a4AC25AB2E6jNn3HqCNJBKnAOQFII0sqLrMC1J1tqKpEqfjwWjJgjUQZCWCNJABnAOwKARpJElHddMG2NOGGNeMcb8X36M0f9kjPmzQR4cci6okQ4L11D3CUEa2D3OAUgKQRpZ0zZMG2O+aIz5uaQ5SUclvSPpLbkruf/aN6q/M9jDRO4EQTqokd7acveXLvU1UBOkgd3hHIAkEaSRRZ2Uebwt6W+stX8St9EYU5ZrZI/388CQc8ZIhUJjjXRQQ10o9KXUgyAN9AXnACSCII2sahmm/XSxL7e6uMRaO2+MWTbGvGKt/VHfjxD5NT3dOGpHEKgJ0kAqcA5Ao402j1uJixv1dVtbW1paWtLq6irtNjKnXZmH6eQqbX8xynJ/DgkjJRqcCdJAmnAOgFxo3ogsd3tbjTyuvxZBGlnXrszjXhevxScAQ0eQBvqKc8DI24jcN1vXjbHt525t7dHS0k2trq7RbiOzehoar4n+j2cGdIEgDQwV54Dciobn3YZpxwXp97S6uk67jUxrV+bRzW93R/saYyaNMWVjzLk2+7XcDoQRpIGB6Ps5AFnSrDd69zdX2nGbII1caNczfdoYc7DD15qU9B9b7eBnzQouWCkaY0q+1i66X1lcGY4OEaSBgenrOQBZFA3U0Zrn7m1tWS0tLWt19Um/pxUAhqJdmD4oN65oJ8Y72OeUpKt+uSqpLGbRwi4QpIGB6vc5AJkU7VmOruuEixtbW3u1tHRPq6sbBGnkRrswPWut/ctOXsgY82oHuxXUeMX3jh4P31s9b4w53cn7YnQRpIGB6/c5ALkQF67bcz3SD7W6ukmQRq60DNOdNqLd7tsGvRtoiyANDN6QzgFIjU6Gu+uMC9JrWl21BGnkTtvpxPuspnpYLigy7FLQK93qBYwxU8aYBWPMwt27dwdykEg3gjQADEMnFxfutLUlVaubBGnkVtJh+oqkol8uSpqXJGNMIVjnR/uY9Mul6AtYa2ettRPW2olDhw4lccxIEYI0ACQh7sLD8HKzEB1eHtPm5l5Vq9LamgjSyK1Ew3QwcocfraMWGsnjLb99zlo7J9d7XUjy2JB+BGkASFonQ+KF9wuMaXPTamlpiyCN3OvnpC0dsdbOxqw7FrPPjv0wIqxtnFbcWq08ekSQBoChiJusJdojHcxq6GJFPUhT2oH8S7rMA2htelo6e7bejWGtVmZmdOvXvyZIA0Ciup20xdnclJaWNgnSGBmJ90wDTVkr1WrSpUvu8cWLLkifOCE7xq8qACSv2aQtG5I2I/uO+SBNjTRGCwkF6WGMdPGiW750SSv/9E+69bOfyT799HCPCwBGWiejd4xpc3NDS0ufJ0hj5FDmgXTxgXplYoIgDQCptSk3zrQL1Jubm1paepogjZFEmEa6BDXSBGkASIH2Y0u7IL1fa2t7ZK1p/lJAThGmkR5BkD5xgiANAEMXHe4uXCsdBOktLS09r7W1vQRpjCxqppEaK48e6da3v83FhgCQSo1t8+bm51St/rbW18cI0hhppBakwvaELARpAEiJsZjbPknS5qZRtfqbPkjzJTdGG8kFQ8fMhgCQVjsD9ebmPlWrX9D6OjXSgESYxpARpAEgrcKzGrpe6c1Nq2p1r9bXRZAGPMI0hoYgDQBpVw/ULkhvaW2NEA2EUeiEoSBIA0CaBb3Rbnlzc4+q1RWtrdFmA1GEaSSOIA0AWTHmLzZc1tpadPpwABJhGgkjSANAdrgg/ZHW1qJjTgMIEKaRGII0AGSHC9K3tLb2ZNiHAqQaYRqJIEgDQHZsbm6qWq1qbW192IcCpB5hGgNHkAaA7KgH6bVhHwqQCQyNh4EiSANAdmxubmpxcVHr6/RIA50iTGNgCNIAMEhxFwX2flonSAO9IUxjIAjSADAI7UbVCLZ3d3onSAO9I0yj7wjSADAI0SAdfTwWs639aX5jY0PVapUgDfSICxDRVwRpABiEVkF6I3Tf3XjQBGlg9wjT6BuCNAAkoVmQbrU95lUI0kBfEKbRXjQcx4RlgjQAJCGuF7pdoI55FYI00DeEabQ2PS2dPVsP0Na6x9PT27sQpAFgkOLCcrNt7cs8CNJAfxGm0Zy1Uq0mXbpUD9Rnz7rHtZpkLUEaABLXLEC3uydIA4PAaB6jylrJmOaPJff44kW3fOmSu0nSmTPSxYtaefSIIA0AiWoVmMdC9zHPJEgDA0HP9CjqoHRjWzhQBwjSAJApBGlgcAjTo6aD0o0d+58927BqZWaGIA0AGUGQBgaLMo9R06Z0Y0fpRxC0g9KOmRndOnGCIA0AGbCxsaHFxVt68uTJsA8FyC16pkdRk9KN2JrpQqGxRvpP/1T26acTO1QAQG82NjYJ0kACCNOjKKZ0o6GGOmx6mhppAMgYF6TfI0gDCSBMj5po6cbWlrsP11BHEKQBIDvqQbq7qcUB9Iaa6VETKd1oKPkoFHaUejCONABkhwvSd/TkyeawDwUYGYTpUTQ93TiudBCoCdIAkFkEaWA4CNOjKu5iwxCCNABkhwvSywRpYAgI09iBIA0A2eGC9D09ebI17EMBRhIXIKIBQRoA0iQ6bXhk64YhSANDRs80thGkASAt2o3EMaYnT9ZUrd4hSANDRpiGJII0AAxfXICOD9VPnmz4IE2NNDBshOlRER69I/KYIA0Aw9R5iJaCIP0hQRpICWqmR8H0dOOELMHELdPTBGkASJUNtQ7SUrX6EUEaSBHCdN5ZK9VqjTMc+hkQV55+miANAEO10WK58fbkyaoWF28zsyGQMomXeRhjJiXVJJWstTMx26f84lFr7fkkjy2XwjMcXrrkbpJWXn1Vt/70TwnSAJAKcaF6Q8Fp+smTDS0ufqiNDXqkgbRJtGfaGFOSJGvtvKRa8Di0vSxp3lo7K6noH2O3woFa0srEBEEaAIYubti76LoNgjSQckmXeZyS65WWpKqkaFguhtZV/WPsVlDaIR+kf/YzgjQApE5ckN4kSAMpl3SZR0HScujxwfBG3yMdKEm6ksAxZVuLUTq2Hwc10pR2AEBKNQvS97SxwTjSQJql8gJEX/5x1Vpbidk2ZYxZMMYs3L17dwhHlyItRunYZoxUKBCkASDVGvu2CNJAdiQdpmuSxv1yQdK9JvuV4y5OlFzvtbV2wlo7cejQof4fYVa0GKVDtVo9YEta+eEPCdIAMFD9GmFjjCANZEzSZR5XJE345aKkeUkyxhSstTW/PBUEaWNM2V+siKgmo3TozBm3nglZAGDAogE6+rjTU+yYgpE73PB3dwjSQIYk2jMdlG34UTpqoTKOt0LrLxhjFo0xnyZ5bJkUGaVDEkEaAAYubmKVZutaiZZ2MGoHkEWJjzMducgwWHfM389LOpD0MWVWaJSObWfPShcvauXRI4I0APRdu8AcPB4LPW5/ql1ft6pW3yNIAxmUygsQ0YFwjfSZM9LWlru/dEkrMzMEaQAYuLgZC5tta2ZM6+tPVK3eJEgDGZV4zzT6xI/S0VAjffGiVn7zN3Xr298mSANA37WaXCW8vvNe6fX1dVWrtwjSQIYRprNserphXOmVR48YtQMABq7VzIX1iwnr66PrHBekq9rY6NdIIACGgTKPrONiQwAYsk56rBsRpIH8IEznAEEaAAZt5wyFO++bbWu0vr6uxcVFgjSQE4TpjCNIA8CwdBqk68tBkN7cpEYayAvCdIYRpAFgGNqN2tEsSH9GkAZyiAsQM4ogDQBJazYUXnRYvJ2jebge6VsEaSCH6JnOIII0AAxLp7MchnuknxCkgRyjZzpjCNIAMGydTCUe9Eg/0eLih9rc3EriwAAMAWE6QwjSAJAWG4oP0XXr66taXFwmSAM5R5jOCII0AKTVzlC9vr6pxcWaNjdps4G8o2Y6AwjSAJB24RppgjQwSuiZTjmCNACkxUaLm9u+vr6lxcV1ca0hMDromU4xgjQApEW72Qo3tLZmCdLACCJMpxRBGgCyY21tr6rVTYI0MIIo80ghgjQApE3z0+Xa2hZBGhhh9EynDEEaALJjbc2oWv2MIA2MMHqmU4QgDQBpVz9trq1tqFpl1A5g1NEznRIEaQBIu3CQtqpWlwnSAAjTaUCQBoC0CwfpLVWrTBEOwCFMDxlBGgCywwXp2wRpANsI00NEkAaA7Fhb21S1epMgDaABYXpICNIAkB1ra2uqVqsEaQA7EKaHgCANANlRD9KMfwdgJ8J0wgjSAJAdBGkA7TDOdIII0gCQHWtra1pcXNTWFqUdAJqjZzohBGkAyA6CNIBOEaYTQJAGgOwgSAPoBmG6X6JB2T8mSANAdqyurhKkAXSFmul+mJ6WajXp4kXJGBekz57Vyu/9nm79639NkAaADFhdXVW1WiVIA+gKPdO7Za0L0pcuSWfP1oP0P/2Tbn3zmwRpAMgAgjSAXtEzvVvGuB5pyQXqS5e0MjGhW7Ozsk89NdxjAwC0RZAGsBv0TPdDKFCvTEzo1s9+RpAGgAwgSAPYLcJ0PwSlHUGQfvrpYR8RAKANgjSAfqDMY7fCNdKUdgBAJhCkAfQLYXq3jHGjdvy7f0eQBoAMIEgD6CfC9C6trKww/B0AZARBGkC/UTO9C0zIAgDZQZAGMAiE6R4RpAEgOwjSAAaFMN0DgjQAZAdBGsAgJV4zbYyZlFSTVLLWznS7fdgI0gCQHZ999pmWlpYI0gAGJtGeaWNMSZKstfOSasHjTrcPG0EaALKDIA0gCUmXeZyS63WWpKqkcpfbh4YgDQDZQZAGkJSkw3RB0nLo8cEutw8FQRoAsoMgDSBJXIDYBkEaALKDIA0gaUmH6Zqkcb9ckHSvy+0yxkwZYxaMMQt3794dyEEGCNIAkB0EaQDDkHSYviKp6JeLkuYlyRhTaLU9zFo7a62dsNZOHDp0aGAHSpAGgOwgSAMYlkTDtLW2IknGmLKkWvBY0ltttieKIA0A2UGQBjBMiY8zba2djVl3rNX2JBGkASA7CNIAho0LEEMI0gCQHQRpAGlAmPYI0gCQHQRpAGlBmBZBGgCyhCANIE0Sr5lOG4I0AGTH48ePdfPmTYI0gNQY6Z5pgjQAZAdBGkAajWyYJkgDQHYQpAGk1UiGaYI0AGQHQRpAmo1cmCZIA0B2EKQBpN1IhWmCNABkB0EaQBaMzGgeBGkAyI7Hjx9raWmJNhtA6o1EzzRBGgCygyANIEtyH6YJ0gCQHQRpAFmT6zBNkAaA7CBIA8ii3IZpgjQAZAdBGkBW5fICRII0AGTHo0ePdPPmTdpsAJmUu55pgjQAZAdBGkDW5SpME6QBIDsI0gDyIDdhmiANANlBkAaQF7kI0wRpAMgOgjSAPMl8mCZIA0B2EKQB5E2mw/TW1hZBGgAyYmtriyANIHdMlhs1Y8xdSbcG/DbPS/pkwO+RFH6W9MnLzyHxs3Trd6y1hwb8HqmSUJsNAIPQtM3OdJhOgjFmwVo7Mezj6Ad+lvTJy88h8bMAAEZTpss8AAAAgGHK5QyIAAAAkmSMKUk6JemetXZm2MeD/CFMtzc77APoI36W9MnLzyHxswDokDHmgqQpSVVJV0KbDkoqSDpvra1F9i9Ya0938NpFa201eGytrRhjJiSdlESYRt9RM41UM8ZMSqpJKrXqUTDGnKPHAb0yxpSstZUm2zr6HQTQHWPMm5Kq1trzkfVTcmH6aGhdSXLBuIPXnbLWzkbWFSVdtta+1JeDB0KomQ4xxkwaY8rGmHO9bE+TDn6WKX+7kPSxdSrUeM5LqgWPY/YrSzqe5LH1ooP/k5LfZzLpY+tWF5+VqaSPrVv+9+f1Jts6+h0E0LN70RU+CBf9ZzNYV+kkSHsEZiSKMO21O2lm6aTawc9SljQf12ClzCm5HkHJfRWY1uNsq8Pfn9PW2jm5/5Ms/36V5Hqb5iVV0/yzSNs/x3KTzbn5HQSywhhT8IvNPpdNn+d7uwvt9gX6iTBd1+6kmaWTartjLYbWVf3jNCqosTE9GN3Bfz0/n9gR9a7l/4nvjV6UJGvtTBc9MMPQyWch+MajmPKfpZ2C2vwOAui71yXNBm2HMaZojLlqjLnqH5eMMdeNMZeDb8B8iC7LfV6Lxphz/lYIv7B/btl/e3Yhsm0ysj3136whHbgAsa6g1ifNdtvTpKAWxxqpJSup8eKPrBkf9gF0qKDWvz/Hpe1e3XLKa3MLav37VTHGVI0xi5IaaiEBIOJ46NvRgtwf69fCbaC1tmqMOS//R7pvY4LHF6y188YYWWvnjDFVuT/i49rQoqRaKKSfMsaU/fOnJC2EthUk/WggPzFyh57pEeaD29UU9xzWVA/LBUVq6zLUK92pe6GGPPV10834k1BN0mVJr/sLf7Kqpha/gwB27Zq1dt7f5qy1J+XKw65H9qtFHi/Lje5RlXZ0EjUzHh7lQzu/mT0VLPiRRLLc0YQEEabramp90my3PU1q6uxY094DekX1hq4oaV5qqKcrhi7YS3Wdsdr/nyzKNezy92m+oLKm1j/LlKRX/O/WSUmZ+8Mg9DsW+zsIYHD8tSPLxpjLbXatttne8f6ha4isLymZSnFHE1KGMF3XLrhl6aTa7mcJhg6a8cuprP8O9dKWFfpqTtJbfvucb3THlf4LTtr9n8xHtl9L8uC61Pb3KxBcpJjUgfXC/zE2Efk2IPgda/Y7CGCwKmp/bVKt3Yt0+s2YMabge8UPyH2rdrKDMA9IIkxv6yC4Zeak2u5n8esvGGMWjTGfDukwO2KtnfVf/82G1h2L2edolv9P/FePtSDQ+T8SUqmDn2VG0lRwAU+HX78Ojf+j7ED43zz8Oxb3Owhg4IpygXq3Ov3Gckpy5R2+TXhJ6b04HynDBYghcSfL6Ek12SPqXaufxfcWHkj8oEbcqPx++eU0lw8BSI/YUZrkysOO7dy9QSFmXbQOutZi3/D6g8aYyUhHRmo7aJAuhGkAAJAoPyxdWe5iw/DkT0flwvCx8NB4ciN3TASjbvjHwcRRs8HU49bamh8y75zcePdzPpz/SFLJ+Nly/beAk3LfCFblrlupRUYWocwDHWE6cQAAAKBH1EwDAAAAPaLMA+iS/8rwlNy40NQGAwAwwijzQC74+rspuYtPwgPtH5SrfTsf1NSF9i9Ya0938NrFyED/8nV7J/0V3wAAYETRM41csNae9xepVKO9xT74Xpe7sCXQzcxWZUnR0Svm5SYkAQAAI4yaaeTNjtkeQzNblUPrKl2MS03vMwAAiEWYRu6FZuZb7vZ5xpg3lf7ZFQEAwJBQ5oFR8LrcOKThMUsvS5K19iV/QeHrcmOXvik3xulLcqUgy/Jjk/rXmo3UXpdUn878uLX2fGjbpFwNd7B9PEsTswAAgPYI08ib45FB909Juhauo7bWVo0x5+UG/Ze1thJ6fMFaO2+MkR/svyqp2GTUjqJC02kbY04ZY8r++VOSFkLbCnKTBgAAgBwhTCNvrvnp0gNzxphJY8z18HTXqk8xG1iWG92jKnU8tfd4ZJSP6DS2p+Sno/WzcnVz0SMAAMgAaqaRe9baOUnLxph2U8NW22zveP/QRY/WGHPVGDPVxQWPAAAgIwjTGBUVuSHuWqm1exFfb92WMaZgrT0p6YBcffbJDsI8AADIGMI0RkVRvuRil0od7jclufIOa+2cn9yloyAOAACygzCNvDkYXeFH3JiU9Eqb5xZi1kXroGst9g2vP+hH8wijzAMAgJzhAkTkgp8evCypGhrGTnKzHhYlHYsMjXdB0kQw6oZ/XPbP3R7+zl84eNmvr/oRPkpyI3OUjDHnrLUzPjhPSqr5EUAW/XJ4ZBHKPAAAyBljrR32MQAAAACZRJkHAAAA0CPCNAAAANAjwjQAAADQI8I0AAAA0CPCNAAAANAjwjQAAADQI8I0AAAA0CPCNAAAANAjwjQAAADQI8I0AAAA0CPCNAAAANCj/x+uebd4mDnNxwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(12,6))\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "uepd.plot_dgm(dgm_true, dgm_tot, m = m, M=M, ax=ax1)\n",
    "uepd.plot_hist(h, ax=ax2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one can see, the expected persistence diagram is much more similar to the true diagram (blue point) than the diagram one would obtain by considering the whole point cloud directly (red crosses)!"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}