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 $\\mathbb{ED}(\\mathbb{X}_n)$:\n",
    "\n",
    "$$\\forall K \\subset \\Omega, \\mathbb{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 $\\mathbb{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 $\\mathbb{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)}) = \\mathbb{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 torus), 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) $\\mathbb{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 allows for 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",
    "from matplotlib.patches import Polygon\n",
    "from matplotlib import rc\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sampling functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_circle(n):\n",
    "    theta = 2 * np.pi * np.random.rand(n)\n",
    "    x, y = np.cos(theta), np.sin(theta)\n",
    "    return np.array([x, y]).T\n",
    "\n",
    "\n",
    "def sample_torus(n, r1, r2):\n",
    "    theta1 = 2 * np.pi * np.random.rand(n)\n",
    "    theta2 = 2 * np.pi * np.random.rand(n)\n",
    "\n",
    "    x = (r1 + r2 * np.cos(theta2)) * np.cos(theta1)\n",
    "    y = (r1 + r2 * np.cos(theta2)) * np.sin(theta1)\n",
    "    z = r2 * np.sin(theta2)\n",
    "\n",
    "    X = np.array([x, y, z]).T\n",
    "\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_noise(N, ndim, scale=1., type=\"uniform\"):\n",
    "    '''\n",
    "    noise sample uniforme on [-scale, scale]^2 or gaussien N(0, scale^2)\n",
    "    '''\n",
    "    if type==\"uniform\":\n",
    "        X_noise = scale * (2 * np.random.rand(N, ndim) - 1)\n",
    "    elif type==\"gaussian\":\n",
    "        X_noise = scale * np.random.randn(N, ndim)\n",
    "    return X_noise\n",
    "\n",
    "\n",
    "def add_noise(X_noise, X):\n",
    "    return np.concatenate([X, X_noise])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute persistence diagram using the Cech filtration\n",
    "def alphacomplex(X, hdim=1):\n",
    "    ac = gd.AlphaComplex(points=X)\n",
    "    st = ac.create_simplex_tree()\n",
    "    pers = st.persistence(min_persistence=0.0001)\n",
    "    h1 = st.persistence_intervals_in_dimension(hdim)\n",
    "    return np.sqrt(np.array(h1))  # we must sqrt to get Cech dgm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expected_dgm(X, n, k, hdim=1, replace=True):\n",
    "    '''\n",
    "    Subsample n points in a point cloud X, k times, and return the k computed persistence diagrams.\n",
    "    '''\n",
    "    N = len(X)  # total nb of points\n",
    "    samples = [np.random.choice(N, n, replace=replace) for _ in range(k)]\n",
    "    Xs = [X[sample] for sample in samples]\n",
    "    diags = [alphacomplex(X, hdim=hdim) for X in Xs]\n",
    "    return diags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Discretization utils to plot expected PD in a nice way.\n",
    "\n",
    "def discretize_dgm(diag, m=0, M=1, res=30, expo=2., filt=None, sigma=1.):\n",
    "    '''\n",
    "    Discretize one diagram, returns a 2D-histogram of size (res x res), the grid representing [m, M]^2.\n",
    "    '''\n",
    "    if diag.size:\n",
    "        h2d = np.histogram2d(diag[:,0], diag[:,1], bins=res, range=[[m, M],[m, M]])[0]\n",
    "        h2d = np.flip(h2d.T, axis=0)\n",
    "        if filt==\"gaussian\":\n",
    "            h2d = gaussian_filter(h2d, sigma=sigma)\n",
    "        return h2d\n",
    "    else:\n",
    "        return np.zeros((res, res))\n",
    "\n",
    "def tohist(diags, m=0, M=1, res=30, expo=2.,filt=None, sigma=1.):\n",
    "    '''Take a list of observed diagrams 'diags' and return an histogram\n",
    "    which is an estimator of E(Dgm(X_n)) for some integer n'''\n",
    "    \n",
    "    output = [discretize_dgm(diag, m=m, M=M, res=res, expo=expo, filt=filt, sigma=sigma)\n",
    "             for diag in diags]\n",
    "    return np.mean(output, axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot3d(X, X_noise):\n",
    "    from mpl_toolkits.mplot3d import Axes3D\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "    ax.scatter(X[:,0], X[:,1], X[:,2], marker='o', color='blue')\n",
    "    ax.scatter(X_noise[:,0], X_noise[:,1], X_noise[:,2], marker='x', color='red')\n",
    "\n",
    "def plot2d(X, X_noise):\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.scatter(X[:,0], X[:,1], marker='x', color='blue')\n",
    "    ax.scatter(X_noise[:,0], X_noise[:,1], marker='x', color='red')\n",
    "    ax.set_aspect(\"equal\")\n",
    "\n",
    "\n",
    "def plot_object(X, X_noise):\n",
    "    if X.ndim == 2:\n",
    "        plot2d(X, X_noise)\n",
    "    elif X.ndim == 3:\n",
    "        plot3d(X, X_noise)\n",
    "    else:\n",
    "        print(\"Unknown number of dimension. Input should be a 2D or a 3D point cloud.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_hist(h, ax=None):\n",
    "    if ax is None:\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111)\n",
    "    ax.imshow(h, cmap='hot_r', interpolation='bilinear')\n",
    "    L = h.shape[0]\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.set_xlabel(\"Births\", fontsize=18)\n",
    "    ax.set_ylabel(\"Deaths\", fontsize=18)\n",
    "    ax.add_patch(Polygon([[L,-1], [L,L], [-1,L]], fill=True, color='lightgrey'))\n",
    "    ax.set_title(\"Expected PD\", fontsize=24)\n",
    "\n",
    "def plot_dgm(dgm_true, dgm_tot, m, M, ax=None):\n",
    "    if ax is None:\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111)\n",
    "    x, y = dgm_true[:,0], dgm_true[:,1]\n",
    "    ax.scatter(x, y, marker='o', color='blue', label='true dgm')\n",
    "    ax.scatter(dgm_tot[:,0], dgm_tot[:,1], marker='x', color='red', label='total dgm')\n",
    "    ax.set_xlim(m,M)\n",
    "    ax.set_ylim(m,M)\n",
    "    ax.add_patch(Polygon([[m, m], [M, m], [M, M]], fill=True, color='lightgrey'))\n",
    "    ax.set_aspect(\"equal\")\n",
    "    ax.set_xlabel(\"Births\", fontsize=18)\n",
    "    ax.set_ylabel(\"Deaths\", fontsize=18)\n",
    "\n",
    "    ax.legend()\n",
    "    ax.set_title(\"Diagrams\", fontsize=24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical illustration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X = sample_torus(nb_points, r1 = 3, r2 = 1)\n",
    "X = sample_circle(nb_points)\n",
    "X_noise = sample_noise(nb_noise, \n",
    "                       ndim = X.shape[1], \n",
    "                       scale = np.std(X),\n",
    "                       type=\"uniform\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_tot = add_noise(X_noise, X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD3CAYAAADhRcqHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAdiklEQVR4nO2d248cR/XHT+0147W9u961HTuJbWxxCeJpPeFiUB6QEf9AAhJ/wOYBpCUixIgQsSAS4QglmDfy+wsi/f6DJA+IJOayRgovvGADCbnYu17P3ryzMztTPFS3p2e2e7qquqq7qvv7kVaz05eZmr6crjp1zvkyzjkBAEAcI0U3AADgLjAQAIBEYCAAAInAQAAAEoGBAAAkAgMBAEhkrOgGpDE/P8/PnTtXdDMAKC03btxY45wfj1vnvIE4d+4craysFN0MAEoLY+w/SeswxAAAJAIDAQBIBAYCAJCIER8EY2yBc/63hHVPEVGDiBY4568kLQN+sbMjXvf3iZpN8X50lGh3l4gxok6HaHZWbFOrifdHjxKNjxfXZqBOZgPBGLtMRFeJ6GLMugUiIs75W4yx8+H7wWVJxgXkz8YG0cgIUbst3rdawgAwJozB1BTRqVP6nz86SvT++0RHjhBxLr5ncpJobk58z8iIeF+rmfk9IBuZDURwo68nrP4uEb0Z/H+LiC4T0VzMMhiInGm1iO7fJ+p2e72As2ftf2+nQ/SlL8ltu7oqtp+cJBobE0ZjdNRu+0A/tqc5Z4goajzmEpYByzQawiAcOtTr+rvO8ZiZ+fl5on/9S/Q0xsaIJibyb1eVcDIOgjG2SESLRERnzpwpuDV+sr8vDAJjYuxfFtbWxPAkyt27omcxNSWMBjCH7VmMBhEdC/6fIaK7Ccv64Jy/zjmvc87rx+MeIyCWZpNoc7PnDJyeLpdxSGJujmhmRvzmQ4eEEdnbK7pV5cCKvWWMzXDOG0T0BhHVg8Xnieit4P+4ZUCDblfMHBw+XHRL3GB3t39ocu+e6Flg9kSPzD2IYMqyHryGvE1EFM5OBDMdDc753+KWZW1DFdncJPr610XXGsYhmdlZ4acYHxfHDKjBXK9JWa/XOXIxBOvrYvYhyzSjTb7/faIf/ED4Pb7whaJbk8zGRjWGXrIwxm5wzutx6+DS8YBWS0z1Fcm//92Lg+h0hKOw0+nFLUxOCsMQwrnwA+zuitdwerLREO/DOIgiEnWnp8XrxgbRQw9hJmQYMBAOs7OT//Dh9m3x2u2KmziMQdAhNBxR5uf733MuDGC7LV5bLWFoTp7U+04VQkNBRLS9LYwf6AcGwkHyNAz37omn+/i4eJoWwcSE+IveoJz3pi/Hxg5ObZomPN57e+hRRIGBcIhm036I8eam6PYfPVqcQZBlLhJCx7kwnK2W+H/OUnjd5GQv1HwEqYzI5nSBMNfBlnEI/QCciyfxiRPuG4c4pqbErMSxY+K33L8v/Aim6XZFz+VrXxP/VxkYiIKxZRjC0GrOy+uIq9VET4hz8XtNT2P+6U/CUFy6ZPZzfQIGoiB2dvq9/iZ47jnxuZwLB1yVMiKnp3szI9vbZkOur18X5yrMcK0S8EHkTKdjPl9gdVWMyU0bHF+ZmhI3c7stehYnTpj53LAX5njokFHQg8iR48fNGodwCDE/D+MQx/i4OOaci/wMUzAmjn0VgIHIgU5HXFQmLtJajWhrS1z0VRpCZGVuThwzU36KqalqGGUYCMucOmWu19BuiycXci/0Cf0UpnoAjIlZqLICH4RFTDxhajXx1EOdA7PUasJQbG1lz8sIe3Jl9E2gB2GB/X0zxiHsMcA42CPsUdy7l/2zyjjkgIEwjInaA9vb4qKFYciPmRlxzLPmgITFfcsCLkGDmHiClLGb6hOffpp9Knp8XAw7yjDTgR6EIbIah2YTxsEVRkfFufjqV/U/Y3e3HLkcRipKMcYuM8aej1m3wBjjjLGbwd/vg+VXg9fFrN/vAlmMwxNPiIux6HoP4CDXr4vehC6c+++XyGQgosI4RNSICuMEHOOcM875BSJ6moTADhHRImPsJgldDG8J4xuy7P+Xv5hrDzDPyEj23oTPRiJrD+K7JKpUE/VEcB4QGI6Q85zz0CA8zTm/MLDeK86ezTZO5bwcXdCqkLU34auRyHqJSongBAVqo8ZgIWlY4gOHDxN98IHevleuwNfgK2FvQhcfjURez7BvBWXwiYiIc/5K0HuYC4xHH4yxRcbYCmNsZXV1NacmynH0aE+4VpVul+jXvzbbHpA/nBP99Kd6+/pmJLIaiFQRnIAHvonAqRmWyL9LQhujD1eFc86cEZF3OpTBYQV6vPSSfjEZn66DrAbiDerd4A9EcBhjM+EGjLHz1PNTEAlfRTjcuEBEXtS0f+EFog8/1NsXQ4pywpj+ufXFSGQyEENEcN4e2HR9YJ/vBL2Imz4I5zz5JNHLL6vvF05hgnJTZiOROZKSc/56zLKLkf9vEdEzafu4ys9/TvTHP6rv9+STRH/4g/n2ADfRHUJm6YXkAUKth8A50S9/qb5fKCgDqgXnQg7xvffU9nPZSOAyHoLOTY74hmrz7rt69TpcHW7gUk5A54S5+hQABhk8yTEnfWuL6LHH1D/aRSMBAxEDjAOIZXmZ6Nlneyebc/F+efnAph98oBcr4VqKPwzEADAOIJZQfOPatZ6RePZZ8b7RiL0IXnpJ+CRU6HTElLorMO741V2v1/nKSj6hEo8+SvTRR2r75Hr4Bl3liL7Kl6hRCFlaInrttaHn4exZ9dD8bje/U8sYu8E5r8eug4EQ6BQJyfMk0vKyeFKFF2N4sc7MxHZxgSUGvdCSF8GhQ6JGhOpX5cEwA4EhRoCqcXjxxRyNg0b3FlggPO5Roj6JIehUl3Kic8g5d/rv4sWL3DbiDMv/HT9uvUkH6XY5X1rqb8jSklgO7BM9/uFxH3wvgeq1Nj5u+XdxzolohSfcf5XvQfzkJ+r73Lljvh2pMCaGF1FSxr7AIIyJ4VzU5/Daa+L9zIz0eVDt7IUSgkVRaR+ETlBTYYdL00EGDGPIUay6i83rDj6IBLw0DktLwjG2tNTvkwD5MHhnaxpn1epUoXBw3jgWlpEfhw6pbV+o1kFS95ZIqXsL3GFkhOjLX5avSRoONbJqrqhSySHG/r7agT59Wj0+wgqIgygdLgw1MMQYQNUKO2EciIx1b4E7qN7wqpGZWamcgVAtGed4BwuUAJXh63vv6Ze608GqcE6w/oBITto+NlFRct7ctNcOAEJGR9WyP0dH7bVlENvCOUQDIjmS+1hhfl5t+yNH7LQDgEFUczVMqJHLYFU4J2BQJEdmH+N0OkR3k2pux4ChBcgbFTmFY8fStzFBHsI5gyI5UmI7plHJtdjettcOAJJQnXrf2LDTjijWnZQ8RSQnDtPCOXt7attPTWX+SgC0UOm5zsykb5MVq8I5CSI5qWI73LBwzkMPyW+LoQUoGpW08EYjfZss2BbOiRPJid3HFs2m/LZ5OX4AGIbKA2121l47iCwL5/AYkZwh+1ihVpPfNo8uGwAyqDgsP/nEXjvyEM6JW5+LcE6rJb9tHg4fAGRRcViePm1vaFzqSMrJSfltVQKoAMgDlZveVo2S0hoIlSIbuordANhG1gl58qSd7y+tgVAZMugoIQGQB9PT8tsaiAg4QGkNhOzsqGqlYQDyZm1NbrsTJ8x/dykNhMp0pcqUEgBFMKcQa2za2V5KAyEbp64ylQRAkcjmET3yiNnvLZ2BUAmMUo19B6AoZH0ROztmyyOWzkB885ty2xVSuh74hYSSd16MjhI9/LDctiaTDUtlILpdouvX5bY1kOIByoyCknde/POfctuZDL8ulYGQrQCF3gMYiqNShypZxipRxMMoVdl7WctZlMYA8ISorMC1az2xIgeEitbW5Cqj3b9v5jovTdn7Vks+tNrxnwxcQVPJ2zayTZC9zitR9l424ElHZRlUkAxK3raRDb82oelZGgMhO+ZSSf8GFcVxqUPZqtYm4nxK44OQCTNFSjeQwnGpQ9ncoW9/m+jPf872XaXwQezsyB20ZlMtBRxUHIelDtfX5UKwO510kWqrPggJ4ZzF4O9qZNkBMZ0syCoNwTgAJRyWOpSd8lSReojDqnBOUFburaCC1PlIVes+MZ2syJSKQ+wDKBOyD7us6QS2hXPOR5bdol6x2kExHW3abbkeRJ5yZQDkgUxgYFZ1uKxOyqEiOAO1JxdIVLQmEmI6REQLnPNXsjRAVvMiLyUiAPJCxn3IuXiIqirah+QyzRkMPd4MK1iniemoCOfITG9ChBeUEdk6qlm0M6wK50S4HPYUEsR0+lARzpHxG6kqaxWOQ1mEwG1sy0TaFs4hxthixDhcpngxHW1khg5e1X1wMIsQuEunk75NFv+bVeGcYPlVxthNxti9yD59Yjq63y8bXu2NgXA0ixC4i0wPWqVk3YHP9zlQanc3/ea/e9czB2XUKIQ4kEUI3EXmshgWJFjaZC2ZWHPv7qloWG8IjAMYgkyMj27iltcGQqb2nnf3lcNZhMaAE9YoMj4GXUe91wZiTCKKwytJPcezCI0AJ2wh6MZBeJ3NmVZX8vOfT09UcQrHswgzE3XCEonfFjWIDiVD+YRMD2JnR+9h6a2Tcns7PYz0zh1Pi9M6nEWYGThhrZB26Obnk6X5SumklJnilJkjdhKHswgzAyesFT74YPj6tTW9+8FbAyFT/yFrogqwQBWcsAUgI6yjU2HKWwMho6AFYV7HqIITtiBkZvR0HJXeOillYtCzRJABC5TdCVsgMiLUOpJ83hqItO7S976H681Jlpf7na6hkcDJchJvhxhpIjmvvppPO4AGZXbCFoRMvhF8EArrASgba2vD1+v4ILw1EACAftIKw+gUjvHWQKQpZEFBC1SNtOLNMsWdB/HWQKSFjXqVgwGAAdLKL+oofmeexQgKvzQooQBt3Pq0fWRIy07zrswcABlJq+4uqx8TxbYuxoH1afvIf3e29QCUjbTAQJ3AQdu6GHHr0/aRIi3oDkF5oGrYGHZnNRBDdTES1qftI0Va5JhMZBkAZcLFHoQVZHQxbBwMAHzGhpPSti5G3PpULQ0ZXYy0uHKduHMApHGwbF5a7pFObpJtXYy49bH7qJKWyo1Ub2ANR8vmpSl56yh9W9XFiFs/ZB8l0opfeFssBriNw9olNnrVmeMgBgR6w2UXU9YfWKZKWh0+qHkDK0RT1K9d65XOc6BsXq2WbX0cTjopZbDhkAFACkfL5qVVWZOpwjaItwYirVq1V9WsgV84WjZvYyPb+ji8vY0QKAUKoWJl87ytKDUxkW09AFo4XDYvTUhKRmjqwD56TSkeJGuBwnC0bN65c8PX6zw0vR1ipImR6oqVAiCFY2XzZB6IlaoolTZl87nP5dMOAFxApt6kTN3KQbw1EDLdJeRjgKog418oTbKWDDLTmIiFAFVBJkqyUj0Imdmkks04AV0cTKzyBW8NhEy9BwRLAVcTq0xz5076NpUKtZbxyG5t2W8HcBiHE6tM8/jjw9d/+KHe53obByHTg3j00VJdA0AVhxOrTLK5mb6NTu+ByOMeBBHR7dvp26AXUXEcTawyicxDULcEo9cGQuZHI6Ky4jiaWGUSGV+bji4nkecGQmZqZ3LSfjuAo1QksWp+Pn0b3R6Etz4IIrmiMLOzqE9ZWRxOrDLF/r5cvI9uCcY8lLUWg38vcM6vBMuucs6vMMYWs1SXknG8dDoiL0MnDh2UAEcTq0who2L/7LP6P9e2stZlInorMALng/dERIuMsZskhHO0kc1Og5BvxXEsscokMr2HX/xC//NtK2udjyy7Rb1q1k9zzi8EhiUTMlM8JRlqAnCAL34xfZssmc1ZhxhDVbIGhg8LJEreExEtMGHFhw1LFomIzpw5M7QBMl2s2VkYCVA+Oh25qf4sBZxzmcUIhh5vhiXuOeevBL2Huciw4wEywjkhst5ZxEOAsiE7dTk9rf8dqT2IiJMxyq3Q70ApKlkBl8OeQuDUJM75/wfbn0/YRwpZ7+zGBsR0QLm4LCF7rSOWEyXVQKTMMrxBRPXg/z5lLc55I/h/MWIcLpPwRYTOyQtE9Hu9pve4d08MI4bx2GMYZoDy0O0S/fWv6dtl9cdaVdYKll9ljN1kjN2L7POdoCdxU1dZK4pstCRmMxRBmrSzyPYMshZvtqqsFQxDDjzbTShrRZmZSd+GCAFTSiwvi4zHMGYgjEqcmSldqrSPyPrepqayfY/XodYhsuHUWZw1lcJ0mjR6IsY5ejR9G5kQgDS8DrWOsrkpd9B2d/VTXyuDyTRp9ESMIzu8MGGHS9GDIJKvHoU6lZKYSJOuUMGWPJFJziKSe2CmURoDITvWkvVXVB4TadKhkQkzKEdGepmVJcqHyJPtbbntTAwviEpkIIjku14ImkrBZJp0BQq25InsjF23a+b7SmUgZK859CJSSEqTXlpST5OuQMGWPJEdXhhzyHPOnf67ePEiV0Fceel/rZbSx1aTbnf4e5n9l5bEAV9ain8PpFlbk7u2t7bUPpeIVnjC/VeaWYyQRkOuhzAxgYdYKlnTpCtQsCVPZHsPJq9rxh2/S+r1Ol9ZWZHevtuVz17b28seaQYkiBZsiXsPUmk00tMJQlRvacbYDc55PW5dqXwQRMJR/pWvyG374x/bbQsIKHHBlryQNQ4y6d8qlM5AEBG9+abcdr/7HYYZwH1UKlKfOGH2u0tpIFTSup97zl47ADDB4cNy262vp2+jSikNBBHR2prcdq++il4EcBeZimkhsoZEhdIaiLm59G1CHnnEXjsAyIJs3tDoqJ3K7aU1EERyisdERJ98Iur7AeASKr0HXeWsNEptIFLKWfYxVrqIEOA7KlnHthTkMhsIxthTjLHLjLHnE9ZfDV4XZfcxyccfy2+rYrEBsInKFLzN69aqcE5An0iO5D7GOHVKflvUiQAuwDnRb34jv71N/VnbwjlEB0VyZPYxykcfyW9rKk0WAF1Uhru7u/baQZTdQAwVzglYGBhOpO7DGFtkjK0wxlZWV1czNpHo9Gn5bVGWDhRJq6WWqq2r2i2LdSclTxHJSdhHWjhHFpWeASKBQVGoDBfyqI5mVTgnQSRHVmzHKKqiOVAEB3kjWzYxJI/r07ZwTpxIzkrcPnmwtydvoZEODvKk3Va73vK6Nq0K5/AYkZwh+1hHNbUbQw2QFyrXZp4lE60K5wxZb1Q4RwXVUgTf+AbRO+/Yaw8AquI2NnIukih1JGUSGxvy2777rrkCoAAMsr+vJgnZbttrSxyVNBCqegGyFaoAUEXF0Tg+nn9KQCUNBJG6JUaUJTCNqo+rCNGnyhqIsTG1sVyzqT5VCkASPhgHogobCCJ1b/D2NhTCQXZUjcPERHExOZU2EETq88kIngJZeOEF9X1k1bRsUHkDQaTefUN8BNCBc6KXX1bfp0hgIEj0CnQ0YYAlBu+Kou8SQ6iGUrswnIWBCNCJdYCRsMDycr92Z6jtubxcZKsyo3qtnD7txvQ6DEQEnQcVjIRBOBcSUlEV8VBlvNHwtiehc42o1DCxCSoxDtDpqFtuxry9dt0iqt157Zr4I+rX9vQMnSa7dC2hBzHAyAjRpUvq+3l47bpJ1EiEwDgUBgxEDO++q7efh9ewe4TDiihRn4Qn6FwLeedZyAADkYDu9Qi18AxEfQ5LS8JzvLTU75PwAB3jUKu5Kb3gYJPcQUelvt2GT0IbxohmZvp9DuFwY2bGiy6abhNVMjrzhPGMV3JQDKZBRAuc81cG1i0Q0Q3qVZV6i3P+DGPsKuf8CmNsMa02RL1e5ysrK5namBXdkw4jocmgZdax1DnT7epPSxZ9nTDGbnDO63HrbOtiHOOcM875BSJ6moiuBsv7tDJcR/cEoiehyaAxcNw4XLrkr3FIw6ouRkQLg4joPOc8NAiDWhnOo3siR0aIXnzRbFuAOzBGdP26+n5TU+4bB6J8dDHC+pNRYzColTG4vVFdDFPontBf/cq+fgHInywdm+1tc+2wSV6zGN8KqlwTUbpWhg1dDFPoGom9Ped7ykCSrC4RH3oOIVZ1MSI88E0kaGV4RZYLhDGRhONCnD1Q54UX1DMyo/hkHIjs62IQYywUywmJ08rwjixGYmyM6JFHiP77X7NtAnYZHc1WwNg340BkWRcjwvrAPn1aGVnaUCRZTvhHHwkDg4rZ7hMOD3XP1dGjfhoHonx0MW4R0TNp+/gK50I7Qzc8e3RUTJPp7g/sktVv1G67GSEpC0KtDfDOOyILVJf33uv5JoAbhBGxWeDcb+NABANhjJGR7N3I8XF1zQ5glnZbnMusOTW+DikGgYEwDOfqpcWibG2JJ9edO+baBOSYmDAj2lwW40AEA2GFTie7hsbJk8JQNJtm2gSS2dkRxzprunWrVS7jQAQDYY3NTTP5/bUaDIUt7t8Xx9aEGC7n5ZREgIGwyNiYuSdKaCjgyMxOoyGOpaqqdhxPPFG+XkMUz32sfsA50e4u0aFD2T9rfFzkddy+DYemKh9/LALUTNHpZPM3+UDJf5471GrmnjTNJtH0tHgKbm6a+cyysr9PtL4ujpUp4zA3l90Z7QvoQeQM52KmwtTTf3pavDYaIujKxHi6DDQaRKdOmffdVC2PBgaiAI4c6T2BTPUqZmZ6/6+uiu+YnDTz2b7QbouwaBsq7Bsb1RzSVaCT5C7drrigTT+Rjh8XfgrGiO7eLbcTrdkkWlsTv3ViwrxxuHJFHL8qGgci9CAKZ2JCdFvbbTsVsefnxesPf0j0s58JX4gJZ2mR7OyIV5vDqampXtBalUEPwhHGx8WTqtm04/z67W+FsZiaEhf91pa40ba2zH+XadbXxazN9nYvbsGWcTh8WBjr8LuqDnoQjjE5KabP9vftBt4MdpnDpyXnYshTq9n77mFsb4vf3WyKtszO5vfdvmde2gCHw1HCICtT8RNpxI3d//GPnj9jb0/0bI4cye4zaTbF7xofF36Ydlsse/TRbJ+rw49+RPTSS9Vz6MpixEAwxhaSCr/E6WYM09IA/YTxE62WiNr7+9/z++7HH49ffuwY0Y0bvQIqnU7vyXvokPD4Hz4s2hwagCNHxHZnzuTT9jTW1sTvwDBiOJkNRFBN6ioRXYxZ90A3gzF2PqqbEV3mc1WpvJiYIHr/fXHD3bkj5viLYn2d6DOfKe77s1DV6UpdMrvDguK16wmr43QzhmppgOGMjBA9/LDoVeztIS1chkZDJGZVebpSF9s+iDjdDCktDZDOxITwEXAunHudTn/AVJW5fVscC4gpZ8PJaU5XhXNc5vBhEXbNuXhaVq1ncekS0b17YjaGc6ITJ2AcTJBVFyONJN2MoVoaQVHb14mEeK/E94AItVq/c7PZFM64snWvt7bEMGt2thqJU0WQVRcjloguRqxuRsIyYIGwjBpRb9p0c7MXMOVLclej0SvKMjHRixHxpf2+YmIW4ykiqjPGngrUsoiELsZFzvnfGGP1Qd2MuGUgH8LeRUjov2i1xF8Y87C1RfTZz+bbtk8/7QVphUFboUEAxcC445k89Xqdr6ysFN2MytJsCoPBuZhi3djoZaHqxEG0Wj0jNToqPmdqChGMRcIYu8E5r8etw2kBQ3nooX5l8mPHkrcNOXnSXntAvsC1AwBIBAYCAJAIDAQAIBEYCABAIs7PYjDGVonoPxKbzhPRmuXm2Aa/wR3K8Dtkf8NZzvnxuBXOGwhZGGMrSVM1voDf4A5l+B0mfgOGGACARGAgAACJlMlAKOeMOAh+gzuU4Xdk/g2l8UEAAMxTih5EtJRdzLqnGGOXGWPP59kmAPLGxn3gvYEIskL/L2Hdg5qYRNQYdgCLJu0EMsauBq9x9TkKQaLNzhtnH497HLbuA+8NhEZNTOeQPIGLjLGbJH5H4aS12Qfj7ONxT8LWfeC9gUjBl/qXMifwac75BclKXnmQ1mYfjLOPx10H7fsA6d5uIHMCF5gQcXBFSyStzT4YZx+Pe644byAs1cT0jojo0LcYY5c9f6J5Q0mOu/Z94LyBsFQTM3dSDN3QExiU9aOgpN9dEr+laNIuOh+Ms4/HXRoT94H3PohoTczI4reJiCI1MAuvf8k5fz3mLzxRb1Dv4ntwAhljocrFLeqd1AtE5EINvrQ2x653DB+Peyy27gMESjlC0MO4RUTnw15TUCvwYmT9erDeibGwZJv71ruGj8c9T2AgAACJeD/EAADYAwYCAJAIDAQAIBEYCABAIjAQAIBEYCAAAInAQAAAEvkfziSz4iX8wwkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_object(X, X_noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgm_tot = alphacomplex(X_tot)\n",
    "dgm_true = alphacomplex(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_points_samples = 100\n",
    "nb_repeat = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resolution parameters\n",
    "m=-0.1\n",
    "M=1.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "diags = expected_dgm(X_tot, nb_points_samples, nb_repeat)\n",
    "h = tohist(diags, m=m, M=M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAF5CAYAAABdm6YkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dX4xc14Hf+d+hWrZH/sNiU7QFz8QjFtdZeHYES8XmbGIEiM0pzbzkZTJNcwE/xLBXzTczFBwyhh/CAPZI5AQSOAG8Ziszb3qg1ZkF8rZmS0YQOBmHzRIHAdY7nnQ1W6IkSiRbRalFkf3v7MM5t+vW7Vt/u+r+q+8HKNSte29XnSK7zv31qfPHWGsFAAAAoH970i4AAAAAkFeEaQAAAGBAhGkAAABgQBNpFwAAACBLjDEVSZOSZK2dT7k4yDjCNAAABWCMOSepKqnid81LaoROKUkq+5skHbPWziVXwlw5IembkuqSDvfyA13+/Uv+viHpjLW23ufPTkpakXSOcJ89htk8MGrDqCSMMSVJV/05s6MtMQDklzHmfUkL1tqn2xwvS7os6aK19nyihRsyY0zFWlsb0XOflnTcWttTmA79XNt/f2NMVdIrkp6L+7fv8rNl/7N1a+2xfsqE0aJlGiNnrT0j9VZJGGM6VRKTav51DwCIt9LpoLW2bow5ISk2bOfMcUkjCdNqbdXvx0q7n7XWzhtjzki6aIypxTQgdfrZuqTDxpirxphXCNTZwQBEJKljJeH/+i8bY16JOd6w1u7LeysKAGRBEboK+EaY3DWwhL5dPTHgU/yhpGljzPSQioRdIkwja6gkACAZd9IuwC4NGkazotL9lJ2stQ1Jc5LODbc4GBRhGplCJQEAo2GMmYnsqseemHHGmJJ/L6fTLssg/BggyY0fGtQluW9yq0MoEnaJPtPIoktyrdPV4KtI/3VeWdLKqAabAEDBHQo/CGbyMMZclBskHszycTioZ/23hEHXu5rct4eTfl9Z0qyki5KCbxP3++c+E1cAP6gv6O53WNIrcV1OfFg+pGbreT0088iMpCP+earhroHRfsR9vN60fz8Nua4jo7zOBH/U7KbRKCjf09pdKMcQEKaRRXGVRFDZz0mKVpYlucopqAT3y42U3tE/25/7A0mL/tyGpAVJU2r2vZtT/IXikKRGaEBl19cNjb4uS/qZpDNy0y1Jfrola+0Jf25QwZbkLhTPxDxXVc3BRZP+foXprQB0Eqo/dgjVQZclTYYbLKy1c8aYuqSnQ9O5NeQHwsl1VaiEx7MYY04bYxblQnm4DrsqV6/VQvsWjTEnwgHXh+N6OJAbY8rGmGlr7VzwWv68laD8Me+519e7qFDd7vdV5OrrofL1/Am1/nsOIrgOlDuehUQQppFFOyoJa+3TcQMTvZlIRV6VdNUYczgmUF+VtF2R+kr0aWvtMR+Oq+ER0/5nqtba836Kv2k1K9iurxt6rsttfuYV/7yLkn4W/JxvTXlVPnD7sp2JXjT8v8mlNv8uAMbXlK/fJuXq0oq6d+s4JmnJGDMTDJLzdVG74FeXVI7+Me/ryxNyLa9BUD8nN5tTtMX3nFyDxSF/3rRcnbsvct4ZNRtUuhrg9aKt9jX/R8SkBlPxZQgEjTVXo681CGttwxgTfl6kiD7TyJxQAI5WEjsqcx9gj4f6oIVHqX8zcu6MXKtL+Cuxi3JdSkp+xpBwRV2XC88/8897Rs1w2/Prhp7rm3It3WGX5Vq365HgPydXGQfPX1X8TCjPxewDgAVr7Qlr7TE/U9IhdZnqzddBx+SmbSv71tl6lxbUdscuqtmdQXL9my/HnLcg1/c3qOvOyde5EVfVY5Ae4PXaPe9iH68XVbfWngndTvjbUNZJCF97kD5appE5fVYSK3KtLpNqvVDUFOkfKBeEoxV/8HhK8f3OGuGQG9ru53W3XyumpXxFUqnNXKMKPX9d0kvGmEuRr2BrvnUCANryc0t37Vvr50GeVe+LurSb0zrocx2eseJITH01qdYwW5YLztFy9RxCQ6/Z6+vtJjSnJWgxz+Ug0qIhTCOLeq4kfLDcJ223FlfkwmdZOyv5oG9z3Gu1uyDElqHP1+34XOphYQAfmuflupHU5S4Gl6218wzIBNCjHS214S4dAWvtCb/I1v5dvFZQ34X79F5qU18FXUqCczsuOtOD4HmSer00BGWPa31HwujmgSzqq5Iwxsz4/s0la+15f2GIC66XtXOwRlmu9bldIG1byfbxukPhR6kHgzKnJV32K2HxdR+ArqLfgAXdOKLn+S5xz0g6vYup14K6tq7moPKO/Y9D3UkGGlTnp8yrJPV6KTsmd+1i8HkGEKaRRT1XEr7SvyjpD7ud7y8kNT+gJnBG7qLRl35edxiCry19S/QJP4Al6E7CnNwABlFVJEz71tpghqATkl7p8gd7u8Aa1Fm1UGiNXb488vx1te8q14tywq+XuNBMUoyZyQjCNDJlgErihKTZmL7I25VlMOWcf+7nJNV9q/KM3AwZg4Thnl93SKaiLUT+gnFMbaa7AoAuTsQMLpwO6kT/bduCpJc6PEe7Vt0TksL9rc+oORd1VLiuPKf4QdyKNIRIzbErgUk1/zjo5/Xa1aFHlM3ZMl6RVOuhPzsSQphG1vRbSVQUP3hkSjtbTKbkW7yttbP+Nmh/435ed1iORXf4C2HXPtcAxkrXOshPmxfdd1k7B2KfkJvxqG3jgJ9eLvz4nCLzNvs6vR593Ui3jCDAz7c5Lxr8X5GrcwPVoE7v8/Wi31gGY2GCQeb9mhzw54KfjQ3wfoaVq5LkZ2hBRjAAEUnqWEnIVYwrHSqJcszPz8t9lReeuzlo2Y5+xVeX+8rymHxf6LiFXUJKat9nup/Xldq/9077FTn2TWPMuXBLkq/wd1wUAYwfH2KrcvVGy8qAXkmuHg1ac4P59l8J/dxxta7+F7TaXvR158XIt3kLcqF1Wq7eOiTpTlw97tcLOO0D7qL8TEXRvtx+3v+Z0HmSa2SJzmc9b4x5zpf/SqTc/b7eadO6WmJdbg7/c34w5jPdvsUM5uSW+3es+D9OXullJpIe/+8ktzAY/aQzxlhr0y4DCi5USQTTFUUrgqCSiFbSwc+X5L5mDFo/5iUdCy1wEoTJYDqlBT/7xWW5CvGibS6NG537VHKV57xcl4+6D/bBiotBeS9Fy9bL6/rnDhZ7Ccp+wr9O+AJWk68kfYV8Qu6CF8zcccVvh1tKSnKtP0OZtxQA+mG6rEAIjAvCNMaGD/VXwqHYB/WyXHj9pqSDXVqrAQAiTAMB+kxjLPivIEsxXxM2/GjzE3KtxrEDXwAAAOLQZxrjIrpSYZwV5XPyfgBIQ0nZnO0CSBQt0xgLvl9xpd0CBMHgGQZ2AEBnflaJy3JjPqrGmMuh1QSBsUOfaYwVP+PGIbVOa3dIkb7UAAAAvSBMAwAAAAOimwcAAAAwoFwPQHz00Uft448/nnYxAKAvW1tbev31129baw+kXZYkUWcDyKvr16/r9u3bJu5YrsP0448/roWFhbSLAQA9W11d1fLysn7/939/Oe2yJI06G0Aeffjhh/rGN77R9jjdPAAgIUGQZqwKAOTDhx9+qDfeeKPjOYRpAEgAQRoA8iUI0t3qbcI0AIwYQRoA8qXXIC3lvM80gGStr6/rxo0bun//ftpFyY2HH36Yfy8AyJF+grREmAbQhxs3buizn/2sHn/8cRkTO6gZIRsbG3rvvfd079497dnDF4EAkHUffPCB3nzzzb6+SaR2B9Cz+/fva//+/QTpHmxubmp9fV2lUonuHQCQA4MEaYkwDaBPRQ7SZ86c0fnz53f9PJubm1pbW5NU7H8vACiKQYO0RJgGkDNzc3Mje+7jx4/v+jnCQRoAkH27CdISYRrACL38svT449KePe7+5Zd393yNRkOXL18eRtFilUqlXf08QRoA8mW3QVpiACKAEXn5ZWlmRrp3zz1eXnaPJelb3xrsORcWFrSwsKC5uTlNT09rbm5Oly5d0vHjx1Wv11WpVHTu3DldvnxZZ86ckSSdO3dOknT+/HlVKhXV63XNBAXxgmO1Wm17X6PR0Ozs7Pb+Uqmkcrmsc+fO6cyZM6rVaiqXyyqVSnrllVf0Z3/2Z3rkkUcGe2MAgMQNI0hLtEwDGJEf/rAZpAP37rn9g6pWqyqXy5qenpYkTU9Pq1araXp6WqdPn1a1Wt0+98SJE9vbZ86cUaVSUbVa1eLiYstzBoG5Wq22/PzCwoJKpZKq1aquXLmimZkZVatVraysqFqtanp6WpcuXVK1WtVTTz2lv/mbvxn8jQEAEjWsIC0RpgGMSLvVV7usytq3SqXS8jiuq0a9Xlej0VCtVtOhQ4dajl29elXlcnnHzwTBen5+frt1W5ImJye3t8vlsjY3N/W5z31uV+8BAJCcYQZpiW4eAEbkS19yXTvi9g9DrVbbEaTD6vX69vaRI0dULpdVqVR2BOfDhw+rXq/v2F+r1bZbwsPCYX1ra4s+0mNno8fH0fu4c6XmZXiiy77osXaPAXQy7CAt8SkEMCI//nFrn2lJeuQRt383yuWy5ubmVK1WNT8/r1qt1hKsjxw5ovn5eTUaDc3Pz6ter+v06dM6f/68VlZWJKmlO8fMzMz2dHi1Wk2XL1/WzMyMKpWKDh8+rHK5rHK5rB/84AdaWFjYfr2f//znqtVqWlpa0tzcnPbt26ejR4/u7s0hozYit/sx+4L99yU9CG2Hb+Gf/VToNhHz+JOh7V5uUnwgb/cYGD93797VjRs3hj73v8nzYgJTU1N2YWEh7WIAY+PXv/61vvKVr/R8/ssvuz7Sb7zhWqR//OPBBx8m7fz585qZmVGpVFK9XtfFixe3u3v0O2vHb37zGz300EMt+5544omr1tqpoRY64/JXZ4dbl+PCcTg8r4ZuDUkf+fvV0P1dvy1Jn/G3vaHtuH3RwP1JtQbv4BaIButegjZQfLsN0t/+9re1sLAQu3AAnygAI/Otb+UnPEdVKhXNz8+rVCqp0Whsz0HN9HdFF9cVIwjN4cD8UeRxEJZv+9ud0PZt6YMt6UNJH/in/LSkz/r7TzwiqSTpUUn7Q9slf8Kn1Bq2P+X3fyZ0LByaJ/z7aNc6vRGzDyiuUbVIB/gUAUCMcFeQAEG66OKCdLgLR9DKHLQwRx/flvSWpJuSbkjvb7ndK/7WkMvYW3I5ea+/L91zt71vS5/5hKTH/C0I1p9RM1iHH2+Gytyu60ec8LG40A0Ux6iDtMQnCAB6QpAuuk4DCze0s/U5SMZBiA6C9HXp/TXplpq3m/5w8HhT0qSkA3J5OdielLR/TTrwhvS5m5J+R83UHbRUf+R/KNq63OkW7gYSvB8CNYoviSAtpfTpMcZUrLW1Nsem5WqnirX2fLIlA4CdCNLjJhqkg5bpoB90Q9K7cgn5plyovinZN6S3/a7g/k1/6lv+/m25MH1ArvH5C5Htx/zTP7YmHaj7LiCPyoX4L6h18KPkLuMPqRma4wJ0XLeOjZjHBGoUR1JBWkrhk2OMqUo6J+lwzLGKJFlr540x5U6hGwCSQJAeB+2mr4sOQAy6dQT9oW/IJeQb0saKC8pvhm7LrfcfrbiH65K++K60vy7pi/722/7+rlzf6qCb9t57UukNaWJVrbOCSK2zgGxG9t1XM1CHz4+GaKB4kgzSUgqLtlhr5+V6j8U5ruZQ57qknZ0WE/Tyy9Ljj0t79rj7l19OszQAhmV+fl5PP/10x3POnDmj559/niA9VjoF6WDGjqBbx025IH1d+nil2eq8LOl/Svo7Sf+vpL+VNv9W+s2K9F8k/Wd/+y+Srt2T7vxPSTVJ1/z5/5+kv1cziL/tX+bjFfda232yb8sl76D7SdByHjcNX/gW9/6A4kg6SEvZWwGxpNagvT+tgrz8spsjd3lZstbdz8wQqIG+RCuzIVRuc3Nzuz6nWq3GrpQYNj09rc3NzY7noKjiFl4JWqY/kguyvkX64w+a/aJ9l2ktyoXp/yG9vyxdkfQrSf9N0mtyYfq/Sfqvcsd+sypt/louTP+dXFPSdTXD9E1/W72nZogP+m0HAyGDQB2e47pTiO70foF8SiNIS9kL05nxwx+2LjYhucc//GE65QFy5+xZ6dSpZoC21j0+e3bgp2w0Grp8+fKuz+lmc3NTn/70p3f1HMiLdisWhu/D80yH5pK2K65LRjBjxy01W6ivSx/dbm2k/h+S/latDdF/J3fOoiQbbY2+peZMIHf9beMDNWcQWVV8eN7UzhbpcKgmPKN40grSUvY6TTXkxjNLrpX6TvQEY8yMpBlJ+tKw1iWO8cYb/e0HEGKt1GhIFy64xy++6IL0hQvSyZPuuImd+76jhYUFLSwsaG5uTtPT05Kk2dlZTU1NaWFhQTMzM7HnBC3VKysrmpmZ6fgazz//vJ544gldu3Zte1+j0dBf/dVf6cknn9S1a9e0d+9eHTx4UC+88IKeffZZXbt2TQcPHtTevXv113/91/rRj37UteUbWRGek7nd/YRa53l+VNJHkrkvTb7d7K58X9I9SWuStqRPPyR9dVl6RG4+jv1yYw03JT0uqSzpdyUdkvTYJ/2Osj/4u5L+gUKzfPjtiUm1TpsXnm86uAUDErutlNhp2XIgPxqNht56661UgrSUkZZpY0xw1bkkV5XI389Hz7XWzlprp6y1UwcOHBhZmdrl9BHmd6A4jHEB+uRJF6D37GkG6RdfHChIS657Rrlc3g7J58+fV7lcVqVSUblc1uzs7I5zarWa6vW6pqendfHixY7P/9Of/lRPPPGEjh49qm984xvb+2u1mvbu3aujR4/q6tWr+u53v6ujR49qZWVFR48e1Z/8yZ9obm5OR48e1VNPPaVajXHT+RYNmcHKg8HqhCW5MPuYNPH55mwcX5QLwV+W9L9K+n3pE1+R/rdHpH8k6Z9I+qf+9k8k/e+Sjkh67AuSvirp9/zPflnNMB3M8PEFSZ+YlJsub7+aU+UF804HQTq6QmK3IB33voH8SDtISymEaT/13ZS/D7wqScHMHX7Gj0aaM3n8+MfSI4+07nvkEbcfQA+CQB22iyAd58qVKyqX3d/f5XI5tntHpVLZXs1wcnJyx/HA5uamrly5ooMHD+44dvToUUnSa6+9ph/96Efb+8PPF/xcp9dA1sUtuR2eai5YLOVRuYT725Ielya+JH1+jwu/j8uF4d+T9ISkiqQp6bGD0j962I2qPyoXqA99UfrEU+64nvTn/56ardP/wN8+v0ea+KLf+dv+tR9Vc9nxYEXEcJDupUU6ug3kSxaCtJTCp8haOydpLrLvcGh7NukyxQmWQP7hD13Xji99yQXpvC6NDCQu6CMddurU0AJ1rVbTkSNHVK/XVS6XVa/XdeTIkR3nLCwsSJJmZmZ07tw5NRqNHV0wgunvnnrqKS0tLe0I1K+//rqOHj26Y//evXt3/T6QBdEp46Lb0TAdnbfZB9jP35D23HMPH5FbLjzolfHbkm5Kn78ptwLiAX/7Qmj7sci+/ZLM59W61Pijal12/DNyQToc+qMLtsQF6ug2kC9ZCdISn6SOvvUtwjMwkCBIh7t2BI+lXQXqcrmsubk5VatVnT59WufPu7WdarWaTp8+veOclZUV1Wo11Wq17Rbqcrm8ve+rX/3q9vR33/3ud/XCCy9Ikq5du6bXXntN3/nOd/TUU0/pa1/7mg4ePKiDBw/q+9//vmq1mq5du6bXX39dv/jFL3Tt2jUtLS1pbm5O+/bt227NRh5FA3UQSoPlvAMPqRmyfd/lR69Ln15p6QmiW3ItzHf89pZcyA76Qke398t36XhMzScqhW5BV5OgdToIzHEt0u2WGI9e/okDyI8sBWlJMlkpyCCmpqZs0OoEYPR+/etf6ytf+UpvJ5896wYhBsE5CNil0q5m9BimXhdkeeGFF/Sd73xHpVJJS0tL+su//MuW7h7d/OY3v9FDDz3Usu+JJ564aq2d6rvQOZb9OrvdtHH31ToNXTCTRvRxMIfdu+5+baX5Yyuhp9iSy8CfVTMfB9vm83Ih+Qv+PmgNDwJ7sB303w62o63P3QYfhhGkkR9pBelvf/vbWlhYiG0F4hMEYDTOnm2dtSPoQz3EPtO70c/Khk8++aR+8YtfaO/evbp7967+9E//dMSlQzqiS20H+8LdPKKPV+UC7n01m6L9HNCfuCN9/ra7qSGtbrmp9CSXix+RNPGImq3MQatzeDt4nej9JyOPw+WPC9bh49H3C+RD1lqkA3ySAIxONDjnMEhLosvG2AlfGsN9pqXmAL/PqBmgwyskhm8fqWVe6s+sSp/5yD/Pp0PPE9zC09yFW5yDWTrCy4dHt9sNnoy+Hy77yKesBmmJTxWAMdNvkMa4C0+PF265bncLL+X9IGbf/cjzhfs6fzLyuFt3jbi+0OrhMZAvWQ7SEp8wAH2y1spkpIW5X2kE6axW/uhHXGAFkIT3339fb7/9dqbr0kws2gIgHz71qU/pzp07ma7U2kkrSDcajdz+8QEAacpDkJb4UxtAH37nd35HN27c0K1bt9IuSl+2tra0ubmZSoVsjCFMA0Cf8hKkJcI0gD48/PDDsSsEZtnq6qqWl5e1Zw9fxGHY4qbSC3B5BQaVpyAt0c0DQIEFQTovFTLyJBqk4x53CtsA4uQtSEuEaQAFRZDG6LQL0nEBmkAN9CqPQVoiTAMoIII0Rmcjsr3RZX/0GIA4eQ3SEmEaQMEQpJGMuLBM1w5gEHkO0hJhGkCBEKQxWnFBuV3LdLt9AMJWVlZyHaQlwjSAgiBIIznRsNxroAYQtrKyonfeeSf39TZhGkDuEaSRnk6BGkA7RQnSEmEaQM4RpJGOTl0+AHRSpCAtEaYB5BhBGgDypWhBWmKJJgA5RZDGeGP1ReRPEYO0RMs0gBwiSGN89TL9Ht1NkD1FDdISf74CyBmCNMZXt9UVJ2KOcZlH+oocpCVapgHkCEEaCLBoDPLhzp07hQ7SEmEaQE4QpDHe2i0A022BGMI10nPnzh3dvHmz8PU2YRpA5hGkgUCnIE2IRnaMS5CWCNMAMo4gDfSyjDmLxiA7xilIS4RpABlGkEb2pDklXbcATcBG+sYtSEuEaQAZRZBGNk10uU9DtzIxoweSMY5BWiJMA8gggjTyIS6sjjLIdnpuAjPSNa5BWiJMA8gYgjSyLy64TqhzgN5N2O30fBPqHOaB0bt9+/bYBmmJMA0gQwjSyLZ2obVdC/EoAm2n5ydAI3m3b9/Wu+++O9b1NmEaQCYQpJE/0WA7yiDdKbAnHegBhyDtEKYBpI4gjfzo1hocDbZx5wyrHN1emyCN0SFIN/FJA5AqgjTyZ0Kt0811upSOYgBiL6/N5R2jQ5BuxacNQGoI0sivuGAbPVbE18a4u3Xrlt577z3q7RA+dQBSQZBGMWRhfmkgGQTpePSZBpA4gjQA5AtBuj3CNIBEEaQBIF8I0p0RpgEkhiANAPlCkO6ODlcAEkGQRnFFBwJyaUUxEKR7wycewMgRpFE8cTNpRI9xiUV+EaR7xycdwEgRpFE80SDdbt5nQjXy6b333tOtW7eot3tEn2kAI0OQRvF0CtLB406t1kC2EaT7l/ify8aYaUkNSRVr7fkOx8vW2tmkywdgOAjSKL6NmO3wgipx20B2EaQHk2jLtDGmIknW2nlJjeBx5HjdH69HjwPIB4I0iikuPPe6DWQbQXpwSXfzOC7X6ixJdUnVmHPO+fuytbaWSKkADA1BGsUXF5g3ItudfgbIFoL07iQdpkuSVkKP94cP+vBcN8YsRs4DkAMEaYyXuPAc3UeIRrYRpHcvUwMQjTEluZbri5JeMsaUUy4SgB4RpDHeCM3IH4L0cCQdphuSJv12SdKdyPEZSc/5gYnHJE1Hn8AYM2OMWTDGLNy6dWukhQXQG4I0MBG5B7KNID08SYfpS5KC1uaypHlpu0W6RTBIMWb/rLV2ylo7deDAgVGWFUAPCNIYP9HgPBG6tTsHyI53332XID1EiX7KrbU1Y8yUMaYqqREaYPiqpMPW2vPGmNPGmLqkSabGA7KNII3xMqH4pcPbLSceF66BdL377ru6ffs29fYQJf7pjgvI1trDoe0dc08DyB6CNMZbOETHXUoJz8gegvRoZGoAIoB8IEhjfEVbm3sJ0gRrpI8gPTp8wgH0hSANRLt2dLqUcplF+gjSo0XLNICeEaSBQLeQ3K7VGkgWQXr0+KQD6AlBGojiEopsI0gng5oAQFcEaQDIl5s3b+rOnTvU2wmgmweAjgjSAJAvBOlkEaYBtEWQBoB8IUgnj24eAGIRpIFR6DQ3NTA4a63effddgnQK+DQD2IEgDQxTdIXE6D4uxdgdgnS6+AQDaEGQBoYpGqTjWqZprcbgrLW6efOmVlZWqLdTwicXwDaCNDBMG222w4+5DGNwBOlsYAAiAEkEaWC42gXpjS7HgN4QpLODMA2AIA2MTC/BmRCN/hCks4UwDYw5gjQwbJ0GHG5Etnv5WaCJIJ09hGlgjBGkgVGKhua41uhOwRpoRZDOJsI0MKYI0kAaurVMA/EI0tnFMGJgDBGkgbRMtNkG2iNIZxst08CYIUgDSZloc99uH7ATQTr7CNPAGCFIA0mIC8jh8NwpXANN1lq98847BOmMI0wDY4IgDaShXWDuZRvjLAjS77//PvV2xhGmgTFAkAaS1iksE6TRGUE6X/jkAgVHkAayILjcsow4OiNI5w+fZqDACNJAmqIBOrwv7jyMO4J0PvEJBgqKIA1kRfhSS8s04hGk84tPM1BABGkgq7jsYieCdL7xqQYKhiANAPlhrdXbb7+tRqNBvZ1TzOYBFAhBGgDygyBdDIRpoCAI0gCQHwTp4iBMAwVAkAaA/CBIFwthGsg5gjQA5AdBungYgAjkGEEaAPLDWqu33npLd+/epd4uEFqmgZwiSANAfhCki4swDeQQQRoA8oMgXWyEaSBnCNIAkB8E6eIjTAM5QpAGgPwgSI8HwjSQEwRpAMgPgvT4IEwDOUCQBoD8IEiPF8I0kHEEaQDID4L0+GGeaSDDCNIA2tuI2cdlPU3WWt24cUMffPAB9fYY4VMHZBRBGkC8uBAdPcblPWkE6fHFpw3IIII0gJ3iQkKmgUkAAB3pSURBVHR430RkP5f4pBCkxxt9poGMIUgD6G5DO8N1dF+nFmwMC0Eaif/ZaoyZltSQVLHWno85XpFUliRr7VzCxQNSRZAGEK9dSI5rmd5os41hI0hDSrhl2gdlWWvnJTWCxxEnfIgutzkOFBJBGkB3nVqeaYlOEkEagaS7eRyXa5WWpLqkavigb7VelCRr7XlrbS3Z4gHpIEgDaK9TX+kNtXbviN5jFAjSCEs6TJckrYQe748cPyJpvzGmYow5nVyxgPQQpAH0pt/ATLAeBYI0orI4APFO0CLtW6pbGGNmjDELxpiFW7duJV86YIgI0gB2j1bppBCkESfpMN2QNOm3S5LuRI4vynX/kL8/En0Ca+2stXbKWjt14MCBkRUUGDWCNID+TOzyHrthrdWbb75JkMYOfYVpY8wzxpjnjDFP+sf/lzHm//H7Hu/hKS7Jz9Th7+f985T8vvnI8Sv9lA/IC4I08mgI1wAMJC4Mh4PyhAjOoxUE6Q8//JB6Gzv0+6lbkfQza+1dY8wzksrW2j+WJGPMP5d0vdMPW2trxpgpY0xVUiM0wPBVSYettXVjTCPo3sHUeCgigjRybFfXAAzDhFpXOWy3aAuhelgI0uim309bw1p7129PS7oYOnY35vwdrLWzMfsOdzoOFAVBGjm362sABtUuRHdqtW53HL0iSKMX/faZDvdxflq+m4bHbxnQAUEaBcA1IDPahWiC9LAQpNGrfj9ph4wxkvR/SJqz1n4gbX+91+j0g8A4I0ijILgGpKpTt464czEogjT60VfLtLX2P8rNsLForf2mMWavMeZfSfoDNQcOAgghSKMouAZkQbT1ud/j6MZaqzfeeIMgjZ71/Ymz1r4U2r4r6c+l7ZYJACEEaRQN14CsIDCPQhCkV1dXqbfRs74/jcaYo9rZAlGSWyr8r4dRKKAICNIoIq4BKCqCNAbVV5g2xjwvV4nWYw6XYvYBY4kgjSLiGoCiIkhjN/ptmb5irf3XcQeMMSywAoggjULjGoDCIUhjt/qdGq/taG0/MAUYawRpFBzXABQKQRrD0Pc80+2WjDXGfH/XpQFyjCCNMcA1AIVBkMawtO3m4ac7qkZ3SyobY/aptc+ckfSUpH839BICOUCQRtFwDUCREaQxTJ36TO+XNKfWFa7a2Scpth8dUHQEaRQU1wAUkrVWy8vL+uijj6i3MRSdwvQla+3rPT7PkjHmuWEUCMgTgjQKjGsACocgjVFoG6bjKlFjzOeC5WP944OSKpLet9a+NpoiAtlEkEaRcQ1A0RCkMSr9DkCcCT+w1i5Za/+jtfY1Vr/COCFIY0xxDUAuEaQxSv2G6U6iK2IBhUSQBmJxDUAmEaQxah0XbfFf4VUlPS1pr9wo7qejp0k6KOniSEoIZAhBGuOEawDyjiCNJHQM09baJUkvSXrJGPMzuVHdszGn1q21d0dQPiAzCNIYN1wDkGcEaSSln+XEn5FUjQ5K8S0X+yRRkaKwCNIA1wDkx9bWlt544w2CNBLRc59pa+3dNsvF3pFkjDH/5/CKBWQHQRrgGoD82NraokUaieqnZVqS5JeSrUoqRQ49Lek/7L5IQHYQpIFWXAOQZUGQvnfvHvU2EtPXbB7GmD+V6y/3v0j6A7mBJ4/67WNDLx2QIoI00IprALKMII209NsyXbbW/pEkGWP+0Fr7anDAGHNUEpP2oxAI0kAsrgHIJII00tTvPNP10Hb0Kz6gEAjSQFtcA5A5BGmkrd8wPSlJfqDJvDHmUuhYdO5RIHcI0kBHXAOQKQRpZEFfYdpa+5Ix5hlJ35RkJb1ijNkyxtyRG9EN5BZBGuiMawCyhCCNrOh7Ng9r7Utyk/hL0pwxZp+kST+5P5BLBGmgN1wDkAVbW1u6fv26Pv74Y+ptpK7fbh4yxjxjjLlijHkutPvgEMsEJIogDfSOawDSRpBG1vQ7Nd6/kht08k25ZWWDifxfM8b88xGUDxgpgjTQO64BSBtBGlnUbzePerACll9CNswMp0hAMgjSQN+4BiA1BGlkVb/dPN4PbUcrzn27LAuQGII0MBCuAUgFQRpZ1m+YPuyXkpXcSG5J25P1U5EiFwjSwMC4BiBxBGlkXV/dPKy1f26M+bn/eq9hjKlLqkiqWWuPj6SEwBARpIHBcQ1A0gjSyINBpsb7I2NMRdJhv+t5a+3rwy0WMHwEaWD3uAYgKQRp5EXfYVqSrLU1SbUhlwUYGYI0MDxcAzBqBGnkSc9h2veJe1ruK71JSQuSfm6t/b9HVLb0WSsZ0/4xcoEgDezeWF4DkAqCNPKm6wBEY8znjDE/lzQn6ZCk1yW9KjeS+8/95P2/O9pipuDsWenUKRegJXd/6pTbj9wgSAO7M7bXAKSCII086qVl+jVJP7XW/lHcQWNMVa6SPTLMgqXKWqnRkC5ccI9ffNEF6QsXpJMnaaEOZLzlniANDMX4XQOQCoI08qpjmPbLxT7TaXCJtXbeGLNijHnOWvuDoZcwDca4AC25AB2E6pMn3f4MBcbUnD3r/uAI/j2ClvtSKROt9wRpYPfG9hqAXdro8bxmBNna2tLS0pLu379PvY3c6dbNw/QyStsPRlkZTpEyIhyoAwRpJ9xyH3SFCVruG41m15iUEKSBoRnfawAGsKHWIL3R5tZ6nCCNvOvWzeNOH89VrE9AEBDDTp0iUEuZbrknSANDNb7XAPQpGpLjtiUXOza2t12Qvq779x9QbyO3+l0BsZPiJMxwS+vJk9LWlrsPt8SOuwy23BOkgVQV5xqAXdgI3Xe6SVtba1paeoMgjdzrFqb7+e3u6VxjzLQxpmqMOd3lvI7HR8oY1/c33NL64ovucalEy7TUvuU+pQqRIA2MxNCvASiyjZj7+JsL0m8RpFEI3bp5nDDG7O/xuaYl/btOJ/hVs4IBK2VjTMX3tYueV1XaI8PPnm2dnSII1ATpnS334dlOpMT/nQjSwMgM9RqAoorr1hHXR9rZ2tqjpaV3df/+Ol/0ohC6hen9cvOK9mKyh3OOS7rst+uSqsryKlrRQEiQdtq13EuJt9wTpIGRGvY1AIXVbuBh+NiEtrbWtbR0V/fvbxCkURjdwvSstfZf9/JExpjnezitpNYR3ztaPHxr9bwx5kQvr4uUZKDlniANjNywrwEotLhZO8J9pK2Wlj4iSKNwOobpXivRfs/tgtaNvEix5Z4gDYxeStcAFEa4j7TV0tI93b+/RZBG4QxzNo9eNNQMyyVFpl0KWqU7PYExZsYYs2CMWbh169aIioksI0gDQH5sbVnV62sEaRRW0mH6kqSy3y5LmpckY0wp2Odn+5j225XoE1hrZ621U9baqQMHDiRSaGQHQRoAsm5i+7a5+ZDq9XU9eGAJ0iisRMN0MHOHn62jEZrJ41V/fM5aOyfXel2KfxaMK4I0AGRZEKKdzc09Wlpa1YMHtEij2LoNQBw6a+1szL7DMefsOA/jiyANAFnWurLh5uaWlpZu68GDTYI0Ci/pbh5A3wjSAJBlEy3bzSDNrB0YD4RpZBpBGgDyIgjSN/XgwRpBGmODMI3MIkgDQF5MaHNzU0tLbxGkMXYI08gkgjQA5IcL0m8SpDGWCNPIHII0AOSHC9JLevDgAfU2xhJhGumIVrj+MUEaAPKDIA2kMDUeoLNnpUZDevFFtwS5tdKpU1r98pe1/PWvUyEDQA5sbm6qXq9rbW2NehtjjZZpJMtaF6QvXJBOnWoG6V/+Ustf+xoVMgDkAEEaaKJlGskyxrVISy5QX7ig1akpLc/Oyj78cLplAwB0RZAGWtEyjeFr0x96WyhQr05NafknPyFIA0AOEKSBnQjTGK6zZ5vdNyR3/y//pdsf2NpyXTuCIP1bv5VGSQEAfQiCNIMNgVZ08yg6a11LcLvHw36toD+05Fqf//E/ln71K+l733PHrZUOH9bqxARdOwAgJ8JBGkArwnSRtZk1Q6VSa0vxsMT0h97h2WddkP7pTwnSAJADBGmgM7p5FFWbWTN04YLbP6qv6MKBOvC970l/8RfSnj1u1o7ZWdlPfnI0rw8AGBqCNNAdLdNF1a6V+OTJZkv1KAShPQaDDQEgPwjSQG9omS6yuFbiJIL0hQsutG9tbbdKM9gQAPKDIA30jjBdZHGtxOGZNobNGNcfO2j99lanplwfaYI0AGQeQRroD908iiraSvzii83H0uhaqM+ebZkxZPUf/kMtf/vbdO0AgNRsxOyLv/xvbm5qcXFRa2troy0SUCCE6aKKthKHu3yUSqPr6hG8tqTV1VUtf/3rzEcKAImLBujw44mY4xMEaWBAhOkii7QSbwfqUQZpb3V1VcvLywRpAEjURg/bYS4GbG4+0OLiMkEaGABhuuiiwZkgDQBjYKPNfZhrod7YMKrXCdLAoAjTGCqCNACkpV2AjgvU4SB9Q2tr60kUECgkwjSGhiANAGmLC9Dh7WZ/6Y2NLdXr72ltrV0XEAC9IExjKAjSAJAV0SC9s2XaBekVgjQwBIRptA5SjHvcBUEaALImLkiHW6RXtba2lVLZgGJh0ZZxd/Zs60IuwfzU/+bftJ7XJigTpAEgC6LhOf7YxsY6QRoYMsL0OLNWajTcQi5BoA4WdvlP/8ktBx6cd+qUC94hBGkAyKKN2JsL0vcI0sCQEabHWTDv9MmTLkDv2ePun3xSunZNevbZ1oDdaGy3UBOkASCLJmJvGxsPqV7f0toadTYwbITpcRdeGTFw9erOgB1aSZEgDQBZEwTnnfs3NvaoXv+YIA2MCGG66KKBN+7xqVOt+559VnrhhdZ9BGkAyIlwi/SWFhfv0rUDGCHCdJG1G1wY9H0Od+E4edL1kf7e99zjw4dbn+vUKYI0AGRes3V6Y8NocfGO1tc3UywPUHxMjVdU4cGFkmtZDgfnYPq7UqnZhePf/lu3P+gz/b3vuZ/91a+0+stfavnv/1724YfTe08AgJ64IP0OQRpIAGG6qMJ9oS9caIbqUN9nSa6VOmhpbjSkf//vpT/4g2aQ/ou/0Opzz2n5j/+YIA0AmRYsES4tLr6p9XUWZAGSQJgusiBQB0Faag3S4fOCY5I7/7//d0nS6vPPa/mf/TO6dgBADrggfZ0gDSSIPtNFFje4MNyHOioys8fq1BRBGgByYmNjQ4uLiwRpIGGE6aKKG1wYTHfXLlCHwvfq1JSWf/ITgjQA5EAzSK+nXRRg7NDNo6iigwvDrc6l0s6uHqHwTdcOAMgPgjSQLsJ0kQWDC4PgHATqaJAOjpVKBGkAyBGCNJA+wnTRtRtsGGP1+99nHmkAyIn19XXV63WCNJAywjQkiQVZACATOg0ebF6yCdJAdhCmQZAGgNTFhegNtV6m3Tnr65YgDWQIYXrMEaQBIG3RIL0Rs+0u1+vrG6rXlwnSQIYQpscYQRoA0tYtSE9sb6+vuwVZNjaYRxrIksTDtDFmWlJDUsVaez7m+IzfPGStPZNo4cYIQRoAsmYjch9sT2h9fUOLizcI0kAGJbpoizGmIknW2nlJjeBx6HhV0ry1dlZS2T/GkBGkASAL4rpz7Ny3vn5fi4tvEqSBjEp6BcTjcq3SklSXFA3L5dC+un+MISJIA0DWtBt8KK2vP9Di4jva2NhMtkgAepZ0N4+SpJXQ4/3hg75FOlCRdCmJQhVOeKGW0GOCNABkWWvrtAvSd7SxsZVekQB0lckBiL77x2VrbS3m2IykGUn60pe+lHTRsu/sWanRaK506JcJX/3yl7X89a8TpAEgB9bXNwnSQE4k3c2jIWnSb5ck3WlzXjVucKLkWq+ttVPW2qkDBw6Mooz5Za0L0hcuSKdONYP0L3+p5a99jSANAJk00XLvZu14lyAN5ETSYfqSmv2gy5LmJckYUwpOMMbMBEGaAYh9Msa1SJ886QL1nj0uSM/Oyj78cNqlA4AC2mhz68XOL4ddkH6HIA3kSKJhOui24UNyI9SN49XQ/nPGmEVjzPtJlq0wgkAtaXVqSss/+QlBGgCGrlto7idUS83p795ksCGQM4n3mY4MMgz2Hfb385L2JV2mQgm6dgRB+rd+K+0SAUDBdFpoRdq5BHinS+2EpA2tra2rXidIA3mUyQGIGFC4jzRdOwBgBOLmho47p/dA7YL0dYI0kFOE6SIxxs3a8S/+BUEaAEaqXaieCO3rfoldW1tTvb5MkAZyLOkBiBih1dVVN/0dQRoARqD94iqtj9uvZhjmgnSdlQ2BnCNMFwQLsgBAUtqF5V7CtkOQBoqDMF0ABGkASFNcC3T7kLy2tqbFxUWCNFAQhOmcI0gDQH4EQXpzkz7SQFEQpnOMIA0AaZqI3Mftax4jSAPFRJjOKYI0AKQlbpaOCcWHa2dtbYsgDRQUYTqHCNIAkIZ2LdBx+5vbtEgDxUaYzhmCNABkRTQ4twvSywRpoMAI0zlCkAaAtEW7cMQ9JkgD44QVEHOCIA0AWTGh+FUPmwjSwPigZToHCNIAkDXRbh1NbrAhQRoYF7RMZxxBGgCyrPUyymBDYPwQpjOMIA0A+fHgwQPV63WCNDBm6OaRUQRpAMgPgjQwvgjTGUSQBoD8IEgD440wnTEEaQDID4I0AMJ0hhCkASA/CNIAJMJ0ZhCkASA/CNIAAoTpDCBIA0B+EKQBhBGmU0aQBoD8IEgDiCJMp4ggDQD5QZAGEIcwnRKCNADkB0EaQDuE6RQQpAEgPwjSADohTCeMIA0A+UGQBtDNRNoFGCcEaQDIjwcPHmhxcVFbW1tpFwVAhtEynRCCNADkB0EaQK8I0wkgSANAfhCkAfSDbh4jRpAGgPy4f/++6vU6QRpAz2iZHiGCNADkB0EawCAI0yNCkAaA/CBIAxgUYXoECNIAkB8EaQC7QZgeMoI0AOQHQRrAbhGmh4ggDQD5QZAGMAyE6SEhSANAfhCkAQwLYXoICNIAkB8EaQDDRJjeJYI0AOQHQRrAsBGmd4EgDQD5QZAGMAqE6QERpAEgPwjSAEaFMD0AgjQA5AdBGsAoTST9gsaYaUkNSRVr7fl+j6eNIA0A+fHxxx9raWmJIA1gZBJtmTbGVCTJWjsvqRE87vV42gjSAJAfBGkASUi6m8dxuVZnSapLqvZ5PDUEaQDID4I0gKQkHaZLklZCj/f3eTwVBGkAyA+CNIAkMQCxC4I0AOQHQRpA0pIO0w1Jk367JOlOn8dljJkxxiwYYxZu3bo1soJKBGkAyBOCNIA0JB2mL0kq++2ypHlJMsaUOh0Ps9bOWmunrLVTBw4cGFlBCdIAkB8EaQBpSTRMW2trkmSMqUpqBI8lvdrleKII0gCQHwRpAGlKfJ5pa+1szL7DnY4niSANAPlBkAaQNgYghhCkASA/CNIAsoAw7RGkASA/CNIAsoIwLYI0AOQJQRpAliTeZzprCNIAkB/37t3T9evXCdIAMmOsW6YJ0gCQHwRpAFk0tmGaIA0A+UGQBpBVYxmmCdIAkB8EaQBZNnZhmiANAPlBkAaQdWMVpgnSAJAfBGkAeTA2s3kQpAEgP+7du6elpSXqbACZNxYt0wRpAMgPgjSAPCl8mCZIA0B+EKQB5E2hwzRBGgDygyANII8KG6YJ0gCQHwRpAHlVyAGIBGkAyI+PPvpI169fp84GkEuFa5kmSANAfhCkAeRdocI0QRoA8oMgDaAIChOmCdIAkB8EaQBFUYgwTZAGgPwgSAMoktyHaYI0AOQHQRpA0eQ6TG9tbRGkASAntra2CNIACsfkuVIzxtyStDzil3lU0u0Rv0ZSeC/ZU5T3IfFe+vW71toDI36NTEmozgaAUWhbZ+c6TCfBGLNgrZ1KuxzDwHvJnqK8D4n3AgAYT7nu5gEAAACkqZArIAIAAEiSMaYi6bikO9ba82mXB8VDmO5uNu0CDBHvJXuK8j4k3guAHhljzkmakVSXdCl0aL+kkqQz1tpG5PyStfZED89dttbWg8fW2poxZkrSMUmEaQwdfaaRacaYaUkNSZVOLQrGmNO0OGBQxpiKtbbW5lhPv4MA+mOMeUVS3Vp7JrJ/Ri5MHwrtq0guGPfwvDPW2tnIvrKki9bap4dSeCCEPtMhxphpY0zVGHN6kONZ0sN7mfG3c0mXrVehynNeUiN4HHNeVdKRJMs2iB7+Tyr+nOmky9avPj4rM0mXrV/+9+elNsd6+h0EMLA70R0+CJf9ZzPYV+slSHsEZiSKMO11u2jm6aLaw3upSpqPq7Ay5rhci6DkvgrMajm76vH354S1dk7u/yTPv18VudameUn1LL8Xaft9rLQ5XJjfQSAvjDElv9nuc9n253xrd6nrycAQEaabul0083RR7VbWcmhf3T/OopJaK9P90RP81/PzyRVpYB3/T3xr9KIkWWvP99ECk4ZePgvBNx7ljL+Xbrr+DgIYupckzQZ1hzGmbIy5bIy57B9XjDFXjTEXg2/AfIiuyn1ey8aY0/7WEqz9z1b9t2fnIsemI8cz/80asoEBiE3dLpp5uqh2LGukL1lFrYM/8mYy7QL0qNvvzxFpu1W3mvG+ud1+v2rGmLoxZlFSS19IAIg4Evp2tCT3x/qVcB1ora0bY87I/5Hu65jg8Tlr7bwxRtbaOWNMXe6P+Lg6tCypEQrpx40xVf/zM5IWQsdKkn4woveMgqFleoz54HY5wy2HDTXDckmRvnU5apXu1Z1QRZ75ftPt+ItQQ9JFSS/5gT951fF3EMCuXbHWzvvbnLX2mFz3sKuR8xqRxytys3vUpR2NRO1Mhmf50M5vZo8HG34mkTw3NCFBhOmmbhfNPF1Uey1r1ltAL6lZ0ZUlzUst/enKoQF7me5nrO7/J4tyFbv8fZYHVHZ7LzOSnvO/W8ck5e4Pg9DvWOzvIIDR8WNHVowxF7ucWu9yvOfzQ2OIrO9SMpPhhiZkDGG6qVtwy9NFtdt7CaYOOu+3M9n/O9RKW1XoqzlJr/rjc77SnVT2B5x0+z+Zjxy/kmjp+tP19ysQDFJMrmj983+MTUW+DQh+x9r9DgIYrZq6j03qWrf0+s2YMabkW8X3yX2rdqyHMA9IIkxv6yG45eai2u29+P3njDGLxpj3UypmT6y1s/7rv9nQvsMx5xzK8/+J/+qxEQQ6/0dCJvXwXs5LmgkG8PT49Wtq/B9l+8L/5uHfsbjfQQAjV5YL1LvV6zeWM5Lr3uHrhKeV3cH5yBgGIIbEXSyjF9VkSzS4Tu/FtxbuS7xQY25cfr/8dpa7DwHIjthZmuS6hx3eeXqLuG8ko/2gGx3ODe/fb4yZjjRkZLaBBtlCmAYAAIny09JV5QYbhhd/OiQXhg+Hp8aTm7ljKph1wz8OFo6aDZYet9Y2/JR5p+Xmu5/z4fwHkirGr5brvwWclvtGsC43bqURmVmEbh7oCcuJAwAAAAOizzQAAAAwILp5AH3yXxkel5sXmr7BAACMMbp5oBB8/7sZucEn4Yn298v1fTsT9KkLnV+y1p7o4bnLkYn+5fvtHfMjvgEAwJiiZRqFYK094wep1KOtxT74XpUb2BLoZ2WrqqTo7BXzcguSAACAMUafaRTNjtUeQytbVUP7an3MS03rMwAAiEWYRuGFVuZb6ffnjDGvKPurKwIAgJTQzQPj4CW5eUjDc5ZelCRr7dN+QOFLcnOXviI3x+nTcl1BVuTnJvXPNRvpe11RcznzI9baM6Fj03J9uIPjk3lamAUAAHRHmEbRHIlMun9c0pVwP2prbd0Yc0Zu0n9Za2uhx+estfPGGPnJ/uuSym1m7SgrtJy2Mea4Mabqf35G0kLoWElu0QAAAFAghGkUzRW/XHpgzhgzbYy5Gl7uWs0lZgMrcrN71KWel/aejMzyEV3G9rj8crR+Va5+Bj0CAIAcoM80Cs9aOydpxRjTbWnYepfjPZ8fGvRojTGXjTEzfQx4BAAAOUGYxrioyU1x10m0tXoH39+6K2NMyVp7TNI+uf7Zx3oI8wAAIGcI0xgXZfkuF7tU6fG8Gcl177DWzvnFXXoK4gAAID8I0yia/dEdfsaNaUnPdfnZuCnwov2gGx3ODe/f72fzCKObBwAABcMARBSCXx68KqkemsZOcqseliUdjkyNd07SVDDrhn9c9T+7Pf2dHzh40e+v+xk+KnIzc1SMMaetted9cJ6W1PAzgCz67fDMInTzAACgYIy1Nu0yAAAAALlENw8AAABgQIRpAAAAYECEaQAAAGBAhGkAAABgQIRpAAAAYECEaQAAAGBAhGkAAABgQIRpAAAAYECEaQAAAGBAhGkAAABgQIRpAAAAYED/P6tl1TLBuD5oAAAAAElFTkSuQmCC\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",
    "plot_dgm(dgm_true, dgm_tot, m = m, M=M, ax=ax1)\n",
    "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 !"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}