Tuto-GUDHI-persistent-entropy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Persistent Entropy and the Entropy Summary Function in TDA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Manuel Soriano-Trigueros - https://personal.us.es/msoriano4/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rd\n",
    "import gudhi as gd\n",
    "import gudhi.representations\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "rd.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will see how to use the variable \"persistent entropy\" (onwards PE)  and the curve \"entropy summary function\" (onwards ES) to compare barcodes using the Gudhi library. Besides, we will explain the intuition behind them and in which cases they can be useful."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please, note that PE and ES are part of the [representations](https://gudhi.inria.fr/python/latest/representations.html) module of the Gudhi library. If you are interested on a specific tutorial about that module and how to use it to create a scikit-learn pipeline, please refer to [this notebook](https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-representations.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Persistent Entropy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is PE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine you would like to use barcodes to perform a statistical analysis or a machine learning task. A common method is to calculate numerical variables from the barcodes and use them in your analysis. For example, you could count the number of bars appearing in the barcode. Not all variables are suitable for every situation. Usually, counting the number of bars has the advantage of having a clear interpretation in terms of the data, but it can be unstable respect to the input in many applications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PE was first defined in its current form [here](https://link.springer.com/chapter/10.1007/978-3-319-29228-1_11), but a precursor of this definition appears in a previous [paper](https://www.sciencedirect.com/science/article/abs/pii/S0031320314002453)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Shannon entropy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to understand what is measuring PE, we need to introduce the Shannon entropy. Given a finite probability distribution, that is a set of numbers $\\{p_i\\}_{i=1\\ldots,n}$ between $0$ and $1$ which satisfy $\\sum_{i=1\\ldots,n} p_i = 1$, its Shannon entropy is the value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "    S = \\sum_{i=1\\ldots,n} -p_i log(p_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Shannon entropy is symmetric, concave and its possible values range from $0$ to $\\log(n)$. The maximum value $\\log(n)$ is reached if and only if all $p_i$ have the same value, that is, $p_i = 1/n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The shannon index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An intuitive way of looking at the Shannon entropy is as a diversity index. Imagine that $\\{p_i\\}_{i=1\\ldots,n}$ are the proportion of the species living in a ecosystem. In that case, $S$ is measuring the amount of diversity. For example, look at the following two ecosystems"
   ]
  },
  {
   "attachments": {
    "blue_red.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAGQCAYAAAD2lq6fAAABhWlDQ1BJQ0MgcHJvZmlsZQAAKJF9kT1Iw1AUhU9TpSqVDlYQcchQnSyIiuimVShChVArtOpg8tI/aNKQtLg4Cq4FB38Wqw4uzro6uAqC4A+Im5uToouUeF9SaBHjg8v7OO+dw333AUK9xDSrYwzQ9IqZjMfEdGZVDLyiG/1UIczIzDLmJCkBz/V1Dx/f76I8y/ven6tXzVoM8InEs8wwK8QbxFObFYPzPnGYFWSV+Jx41KQGiR+5rrj8xjnvsMAzw2YqOU8cJhbzbay0MSuYGvEkcUTVdMoX0i6rnLc4a6Uqa/bJXxjM6ivLXKcaQhyLWIIEEQqqKKKECqK066RYSNJ5zMM/6PglcinkKoKRYwFlaJAdP/gf/J6tlZsYd5OCMaDzxbY/hoHALtCo2fb3sW03TgD/M3Clt/zlOjD9SXqtpUWOgNA2cHHd0pQ94HIHGHgyZFN2JD+VkMsB72f0TRmg7xboWXPn1jzH6QOQolklboCDQ2AkT9nrHu/uap/bv3ea8/sBgeFyrUe241UAAAAGYktHRAAAAKEA/wz7Y6UAAAAJcEhZcwAALiMAAC4jAXilP3YAAAAHdElNRQflCRUJMgbrSC+GAAAAGXRFWHRDb21tZW50AENyZWF0ZWQgd2l0aCBHSU1QV4EOFwAAIABJREFUeNrs3Xl8VYWd/vHnnHvunp1ACBiQJWETkCWW1gUE6lrFCj/LULVoO44g4zhKSytVYGq1IgydVtFpi8aNRaVSVwQ1UUFREGqUAFFUiBK2kD13v+f3xw0RJgioKNzweb9eeRHvcpDvvclznnPPYti2bQsAAAAAcNIxGQEAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAAAohAAAAACAk43FCAAAQDIzDKPle9u2GQgAUAgJEwAAAAA4EnYZBQAAAAAKIQAAAACAQggAAAAAoBACAAB8U1VVVSouLtZtt90mwzB01VVXMRQAOAEk1UllqqqqVFpaqldffVV33HGHrrzySj366KO8igAAnMDi8biys7MZBABQCAkTAABONoZh6KabbtIFF1yggQMHatasWWpoaGAwAEAhJEwAADgZCuG8efMYBABQCAkTAAAAADhRcFIZAAAAAKAQAgAAAAAohAAAAAAACiEAAAAAgEIIAAAAAKAQAgAAAAAohAAAAEfw+9//XoZhyDAMPfDAA3rsscda/vu2225jQABAISRMAAAAAOC7ZDECAADwbZs+fbqmT5/OIADgBGPYtm0zBskwjJbvGQkAAGQ4AJwMOIYQAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAAAKISMAAAAAgJOTxQgAAMDXsX79ei1ZskTLly9Xfn6+xo4dqzFjxsjn8x3V85uamvTqq69qxYoVWrNmjbp166bvfe97uuyyy9S9e3cGDADfAcO2bZswkQzDaPn+OI8EAIATXnFxsUaOHNnq9ilTpmju3LlyuVyHfX44HNbkyZO1YMGCL10/GDRoEBkOAG25EBImAAAkn/r6eg0bNkx+v19FRUUqKChQQ0ODFixYoKlTp+qll17Seeedd9hl/POf/9SgQYM0Y8YMTZo0SdnZ2QqFQnrnnXd07rnn6pZbbtGcOXPIcAD4lh23Ywjr6+s1ZcoUFRYWauPGjYpEIqqurtacOXN07733qqSk5IjLKCsr04IFCzRjxgzt3LlT0WhUjY2NKi4uliQ9/vjjvMIAABxjpaWlKisr0913362+ffvKsixlZGRo8uTJGjZsmBYtWnTEZUQiEUnS5ZdfrpycHDkcDvl8Po0YMUJjx45VY2MjgwaAtlwICRMAAJJTeXm5JGnAgAEH3e71ejVhwgQVFRWpqanpsMvo1auXCgsL9fTTT2v37t2KxWIKBAJ6/fXXtXTpUl1yySUMGgDaciEkTAAASE6VlZWSpMzMzFb35eXlSUrsCXQ4aWlpWrZsmZxOp3JycmRZlnw+n+bMmaOVK1fqggsuYNAA0JYLIWECAEByamxsVEFBgUyz9WqE2+2W9MVePIeze/dubdq06aDbNm7cqG3btikcDjNoAGjLhZAwAQAgOfn9fpWXlysej7e6LxQKSZKcTudhl/Hhhx9q0KBB6tevnyorKxWJRNTQ0KCioiLdeuutuvfeexk0ALTlQkiYAACQnHJzcyVJ1dXVre6rqKiQJKWmph52GftPAHfDDTeoY8eOsixLfr9fZ599tu6880798pe/POKeQgCAJC6EhAkAAMmpoKBAUuIEcQcKBAJauHChJk6ceMTrCR8q//fbf+mIQCDAsAGgrRZCwgQAgOTUv39/FRQUaNq0aSorK1M0GlVNTY3mz5+vNWvWaPz48UdcRq9evSRJ991330GXjnrjjTc0ffp0FRYWKisri2EDQFsthIQJAADJKS0tTfPnz9fatWvVr18/OZ1OZWZmaurUqZo0aZJGjBjR8thAICDDMDRr1qyDljFy5Ehdcskl+s1vfqPc3Fw5nU6lpKTonHPO0e7du3XPPffIsiyGDQDfMsPe/1HacfDKK69o9OjRrW6fNGmS5s2b13JymUAgIJ/Pp5kzZ2rGjBktj6urq9OVV16pZ5999pDLLykp0fDhw49uEIbR8v1xHAkAAElj3bp1Wrx4sVauXKkePXpo3LhxGjNmjPx+/0GF8FAZLkkNDQ166aWX9Nprr2n16tVKS0vT6NGjNWbMGJ122mlHvzJDhgNAchZCwgQAAHzjlRkyHACStxASJgAAgAwHgOPDZAQAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAD4ZtavX69p06Zp4MCBGjdunBYtWqSmpqavtIxgMKjnn39ekydP1pAhQ/Tv//7vWrlypcLhMAMGgO+AYdu2fbzDZMmSJVq+fLny8/M1duxYjRkzRj6f7yuFySuvvKLnn39eb7/9tn7wgx/o0ksv1fDhw+VyuY5uEIbR8v1xHgkAACe84uJijRw5stXtU6ZM0dy5c48qf+vq6nTzzTdrwYIFre4rKSnR8OHDyXAAaMuFkDABACD51NfXa9iwYfL7/SoqKlJBQYEaGhq0YMECTZ06VS+99JLOO++8Iy7nzjvv1PTp0/Xb3/5W//Zv/6aOHTsqGAxq3bp18vv9KiwsJMMBoK0WQsIEAIDktHr1ap111ll69dVXde6557bcHggENHLkSPXu3VsPPfTQYZdRUVGhLl26aPz48SoqKpLb7f76KzNkOAAkXyEkTAAASE4PPfSQrr32Wu3du1ft2rU76L4///nPuvHGG9XY2HjYwz9efPFFXXTRRa3WA8hwAPhuHbeTypSXl0uSBgwYcNDtXq9XEyZMUFFR0REPTP/ggw8kSdddd903KoMAAODoVVZWSpIyMzNb3ZeXlycpsSfQ4Wzfvl2SlJ+frxUrVmjcuHEaMmSIbrvtNr333nsMGQDaeiEkTAAASE6NjY0qKCiQabZejdi/gTYSiRx2GQ0NDZKkl19+Weeff76WLl2q9evX64477tDpp5+uV155hUEDQFsuhIQJAADJye/3q7y8XPF4vNV9oVBIkuR0Og+7DI/HI0latmyZNm/erFAopEAgoPXr12vAgAG68cYbFQwGGTYAtNVCSJgAAJCccnNzJUnV1dWt7quoqJAkpaamHnYZOTk5kqSpU6eqV69ecrlc8ng8GjRokGbMmKGysrKWvYkAAG2wEBImAAAkp4KCAklSaWnpQbcHAgEtXLhQEydOPOL1hPPz84/49xx4shgAQBsrhIQJAADJqX///iooKNC0adNUVlamaDSqmpoazZ8/X2vWrNH48eOPuIzevXtr5MiRmj17trZs2aJwOKxgMKgNGzZo1qxZGjx4cMvGYwDAt+e4XXairq5OhYWFSk9PP+R1CJcvX67zzz//sMsIhUK66KKL5Pf7dc8996hbt26Kx+PatGmTJk6cKMuy9Oabbx7VGUg5ZTUAAEfvlVde0ejRo1vdPmnSJM2bN68lewOBgHw+n2bOnKkZM2Yc9NgNGzZo8ODBh1z+qlWrdOaZZx7dygwZDgBfm3W8/uK0tDTNnz9fo0ePVr9+/VqFyYgRI1r++8vCxO12a86cORo8eLCeffbZQ4YJl6MAAODYGzVqlNauXavFixdr5cqV6tGjh8aNG6cxY8YcdfYOGjRImzdv1hNPPKEXXnhBknTRRRdp3Lhx6tOnD0MGgO/AcfuEcL9169YdMkz8fv8RC+F+W7Zs+cZhwtZFAACSdGWGDAeA5C2EhAkAACDDAeD4MBkBAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAQCEEAAAAAFAIAQAAAAAUQgAAAAAAhRAAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAAAAIBCCAAAAACgEAIAAAAAKIQAAAAAAAohAAA4JtavX69p06Zp4MCBGjdunBYtWqSmpqavtaxAIKDLLrtMhmHoH//4B8MFgJOlEBImAAAkn+LiYg0ZMkSzZ89WaWmpli5dqgkTJmjatGkKh8NfeXmLFi2S3+9nsABwMhVCwgQAgORTX1+vKVOmqLCwUBs3blQkElF1dbXmzJmje++9VyUlJV9peVu2bFFRUZGmT5/OcAHgZCmEhAkAAMmptLRUZWVluvvuu9W3b19ZlqWMjAxNnjxZw4YN06JFi456WcFgULfffrvmzZunrKwshgsAJ0shJEwAAEhO5eXlkqQBAwYcdLvX69WECRNUVFR01Id/PPXUUzrrrLM0ZMgQBgsAJ1MhPCnCxDCO/RcAAMdZZWWlJCkzM7PVfXl5eZISewIdydatW7V48WJdc801DBUATrZCSJgAAJCcGhsbVVBQINNsvQrhdrslSZFI5LDLCIfDmjVrlv7whz8oJSWFoQLAyVYICRMAAJKT3+9XeXm54vF4q/tCoZAkyel0HnYZy5Yt0znnnKPTTjuNgQLAyVgICRMAAJJTbm6uJKm6urrVfRUVFZKk1NTUL31+Y2OjfvKTn+hf//VfZRhGy9f+5e6/hNSbb77JsAGgrRZCwgQAgORUUFAgKXGCuAMFAgEtXLhQEydOlM/n+9LnR6NRhggAJ3shJEwAAEhO/fv3V0FBgaZNm6aysjJFo1HV1NRo/vz5WrNmjcaPH3/Y56enp8u27VZf+88vsGzZMtm2rR/84AcMGwDaaiEkTAAASE5paWmaP3++1q5dq379+snpdCozM1NTp07VpEmTNGLEiJbHBgIBGYahWbNmMTgAOEFZxzNMRo8erX79+h1036HCxOfzaebMmZoxYwavGAAAx9moUaO0du1aLV68WCtXrlSPHj00btw4jRkzpuXkcAAACiFhAgBAGzV06FANHTr0sI/xer2ybfuoltexY8ejfiwA4NgxbH77JgZxwIXfj8lIvo0LyfNSAQDw7Wc4AJxETEYAAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAABRCAAAAAACFEAAAAABAIQQAAAAAUAgBAAAAAF/OYgQAAOBkZzxy7JdpX81cAZz4+IQQAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAANoQzjKKE4NhHPtl2jZzBQAAAA6DTwgBAAAAgEIIAAAAAKAQAgAAAAAohAAAAAAACiEAAAAAgEIIAAAAAEh2XHYCAAActfXr12vJkiVavny58vPzNXbsWI0ZM0Y+n++Izw2FQiotLVVJSYnWrl2rPXv26IwzztCFF16oM888U06nkwEDwHfMsO3v/mJtJ2KYGAdcB++YjITr6jEvAGhjiouLNXLkyFa3T5kyRXPnzpXL5Trs85988kldccUVh7zvlltu0V133fW1cvxYZLjxyLcQQ1fzngFAITypwoSCw7wAoK2qr6/XsGHD5Pf7VVRUpIKCAjU0NGjBggWaOnWqXnrpJZ133nmHXcaLL76opqYmDR06VNnZ2XI6ndq2bZtmzpyphQsX6u2339YZZ5xBIQSA79B3egxhfX29pkyZosLCQm3cuFGRSETV1dWaM2eO7r33XpWUlBxxGSkpKXrqqaf06aefqqGhQaFQSOXl5ZowYYLmzp2rDRs28KoCAHCMlZaWqqysTHfffbf69u0ry7KUkZGhyZMna9iwYVq0aNERl3HhhRdq7Nix6tq1q/x+v1wul/Lz83XrrbdKknbv3s2gAaAtF0LCBACA5FReXi5JGjBgwEG3e71eTZgwQUVFRWpqavpay3a73ZKkTp06MWgAaMuFkDABACA5VVZWSpIyMzNb3ZeXlycpsSfQV2Xbtp544gldd911rdYPAADfvu/0LKNHGyZHc3IZwgQAgO9OY2OjCgoKZJqttyXv3ygbiUS+cn4/+uijKikp0WOPPSbL4uTnANCmCyFhAgBAcvL7/SovL1c8Hm+V46FQSJK+0kndotGo/vrXv+rFF1/UI488og4dOjDkJMZJeZgh86QQEiYAALRhubm5kqTq6mq1a9fuoPsqKiokSampqUe1rEAgoNmzZ6usrEwPPvigsrOzGTAAHCff6TGEB4bJ//V1wuT3v/+9SkpK9OCDD6pjx468mgAAfEsKCgokJU4Q93/zeOHChZo4ceJRHfJRU1OjW265RZ988okeeOAByiAAnEyFkDABACA59e/fXwUFBZo2bZrKysoUjUZVU1Oj+fPna82aNRo/fvwRl1FZWalrr71WkvSnP/3pkOcUAAC04UJImAAAkJzS0tI0f/58rV27Vv369ZPT6VRmZqamTp2qSZMmacSIES2PDQQCMgxDs2bNOmgZ//jHP/T000/r/vvvV3p6ugzDOOjroYceYtAA8B37To8h3B8mo0ePVr9+/Q6671Bh4vP5NHPmTM2YMaNVmEjS/fff3+rvePDBB3XNNdfwygIAcIyNGjVKa9eu1eLFi7Vy5Ur16NFD48aN05gxY1pODgcAoBASJgAAtFFDhw7V0KFDD/sYr9cr27Zb3X799dfr+uuvZ4gAcAIx7EP9xj4ZB2EYLd8fk5EcsLxjpi2/VMwLAHAcM5xT1H/D14D5MUPmmbRMAQAAAABOShRCAAAAAKAQAgAAAAAohAAAAAAACiEAAAAAgEIIAAAAAKAQAgAAAACSncUIgDaGazoCAADgKPEJIQAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAAEAhBAAAAABQCAEAAAAAFEIAAAAAAIUQAAAAACiEAAAAAAAKIQAAAADg5GAxAgA4AsM49su0beYKAACOOz4hBAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAAAAAoBACAAAAACiEAAAAAAAKIQAAAACAQggAANq89evXa9q0aRo4cKDGjRunRYsWqampicEAwAnIOpHDZMmSJVq+fLny8/M1duxYjRkzRj6fj1cNAIATVHFxsUaOHNny36WlpVq6dKmmTJmiuXPnyuVyMSQAOIGckJ8QFhcXa8iQIZo9e3ZLkEyYMEHTpk1TOBzmVQMA4ARUX1+vKVOmqLCwUBs3blQkElF1dbXmzJmje++9VyUlJQwJACiEhAkAAG1RaWmpysrKdPfdd6tv376yLEsZGRmaPHmyhg0bpkWLFjEkAKAQEiYAALRF5eXlkqQBAwYcdLvX69WECRNUVFTEsYQAQCEkTAAAaIsqKyslSZmZma3uy8vLk5TYEwgAQCEkTAAAaGMaGxtVUFAg02y9euF2uyVJkUiEQQEAhZAwAQCgrfH7/SovL1c8Hm91XygUkiQ5nU4GBQAUQsIEAIC2Jjc3V5JUXV3d6r6KigpJUmpqKoMCAAohYQIAQFtTUFAgKXGCuAMFAgEtXLhQEydO5HrCAEAhJEwAAGiL+vfvr4KCAk2bNk1lZWWKRqOqqanR/PnztWbNGo0fP54hAcAJxjqRw6SoqEgFBQVqaGjQggULtGbNGs2cOZNXDQCAE1BaWprmz5+v0aNHq1+/fgfdN2nSJI0YMYIhAQCFkDABAKCtGjVqlNauXavFixdr5cqV6tGjh8aNG6cxY8a0nBwOAEAhJEwAAGijhg4dqqFDhzIIAKAQEiYAAAAAcKIybNu2GYNkGEbL98dkJAcs75hpyy8V82KWzBTAMchwAGjLvo3qZjHWJAoWAo95MUtmeoKGCQAASE4mIwAAAACAkxOfEAIAgKTGp94AQCEkTAAAAADgK2KXUQAAAACgEAIAAAAATibsMgoAAIAjisVi2r1rpxobG5WSkqoOOTkyzaP/bKGhoUF7du2SDKl9hxylpKQwVIBCCMKEMAEAJJ/dNXtkGIbap2e3+X+rbdt6+83VevDPf9SH77zccnvfsy/WtTf8hwYXFh72kl011dV64vHH9NT99ygSqJMkOb1pGjfpl7rip1cqIzOTNxRwHHFhesKEMAEA4Ct4e8ta/eKF30iGoQUX36Uz8oe26fxe8tij+t8ZN+ria27SqAsvVkZmpvZVVemlZ5Zp5cL79Z9zFujSy8ceel1n1y796vqfq6GmStfeNE35vXpJksq3bNZDf5ytlIx2mv3AAnXIyTkp30uhUEj19fWyLEupKSlyWMf2s5pYNKr6hgZFo1GlpqbK7XbzAwwKIWFCmBAmhAkAfF2vvF+sq17+rQKRmGxJfo9Dj466UyP7DW+T/97169bplitGa+p/P6SLxlwm2zZkS3rztWKtffMNvf3Kcu3+tEz3PbtKffuddnB+xGL67S03afdn2/WH+X9R+w45sm0pbkvvvbtWzy9dopKlRfrejybod3PmyeFwnDTvo7KyMi17+u+q2LZN2e07KBwJa8+uXTpnxLm68MILldOx4zda/q6dO/Xiiy/qzdWrlJPTUYZpaNfOnera9VT9eOxY9enThx9mUAgJE8KEMCFMAOBoba+q0MbPN+lnq25XPCLZcSkekxwuyeEy9MhZv1OfTr3UJSuvTf27fzVlkrxej2bcPVemaSoWs2XL0NJFj+mB26e0PG70Tydr+u/uPOi5mzZu1ORLztSf/l6s/qcPkiR9VvGZOnY6RVvLt+j6i7/X8tgHnn9Tvfr0bfPvo2g0qocfekg7duzQ2cOHy+/3KxaLtXx9/tlnKn7lZY37yXj98Ic/POzeU4di27ZWrlih5S++oAsuulg5OTmKxWKKRqOKxWKqrt6n10teU8/8fE285hpZFkePgUJImBAmhAlhAgBH/L04bulEvVVXrmjIlh1NZHgsbMvpM2Q4DFlu6ZyMflp42V+/8u/dE1VdXa1+XNhLMx9YqLPPHSlJeu7ppVrx7D/0wWvPy7ZjLY91p2Zr2dsfyOPxtNz2zNNL9eCcO/RkyTtyOp2q3leln5w9WKedOUqb3nlDwbo9LY+9ee6DuuTHl7fp91E8Htf8++5Tamqq+p122kHZvT9j93//1upV6j9goH4yfvxX+juWLF6sXbt2afiIEYrH4y3L27/8/d9v/OADuVxOTb5hylc6lwPaJtbkvsMw+c/i3+qtunLJkIz9P3tRW6YjERzXrr2tTYbJuyuf1swHFrb8wnnxmb+3hMmB3njmCQWn335QmHxYvlnpud3Uu/mTw+p9Vbr2gjNbwuRA5Zs3t/lCGI/H9cD99ys1NVUXXnyxYrGYIpHIQWGS1a6dxlw+Vq+XFKt6376vHCZPLFmiXbt26WfXXKt4PK5IJHJQmLjdHp07apQ2fvCBHrh/PmECoM3btm+73gl+KEkyDEO2mrelG4liaFqSLWlVoEzbqyvUNatLm/h3h0NhxSNBpaSmttxWse1TvV/yTKvHhur3KhwOHZThTQ2N6pDXQ06nU5JUW1urSKBWG17+e6vnNzTUt/n30RtvvKFoJKI+ffsqHA4rsGunakvfU2THDsnplHVKnpzduivudmvQkKFa9frrysvL0w/OPPOolr969Spt375NPzzvfEWjUYVqa1X9wfsKfLxV8WBQZvsOchf0kjIy1L1HD73z9hq98cYbGj48ufdQCwQCWr/+XQUDQeV16aLu3bt/Zxuro9GoPv74Y1Vs3y6P16PBg4fI6/VSCEGYECaECWFy8oUJ0Fa9t+t9/WPz81LzdlrDkIzmoxMMUzIto+WqzoakRzct1qX5F2lAh9OS/t/uT0mROzVblZ9/rkFDCyVJXbv3OORjM04pkM/nP+i27Pbt9WnpGjU2NMifkqL09Aw5nB7FIsFWz2/foW2fB6CpqUkyc0c4AAAgAElEQVSPPlykidf+XOFwWPveeVu1ix+Tx+OSx+WUbClc+q7qbEOO8y6S0TNfg4YM0b1/+h/1O+00paenH3b5tbW1euyRR3Ttz3+hcDisqvXvqmrRY3LbUblciXWoyKZS1ax4Xs4LL5XZp6/69O2nvz5wvwoLC+Xz+ZJyrqFQSLfd+htt2/SeLMtSQ0OjMjt11TW/+FededZZcrlc3876bTis1atW6aG//VXVO7YpJcWvaDSqU/sO0h133ZV051ugEBImhAlhQpgQJgC+xNu71+kvVc/IMCTbTgS1YRgyLFvxqJHIczuR7TKk+6qeVsesnDaR4V6vVxdfdZ2efPhvGjH6h/L5/Tpr+LlaOugcfbzh9YMee+mV17bakNa3f39Fgg1aVfKqzv/RpcrMytI1v/69/va7Ww56nOn0qF///m36ffTuunUaWniGYrGYGsq3qOHJhcrKTJPP75PLackwTEUiETUFgqpb/qzCQ4cpNnCQ+vQ7TS+vXKmx48Yddvkvr1yhc0eOUiwW065XXlbwxWeUnZ4qny9DltOSbFvhcESNTQHVvvScIk6XYp07q8up3fTuunU6+5xzknKuxcXF2r65VFkZiXWctBS/7FiT7v3Df2lRXnfd/MtfqXfv3sf079y8aZP+e8492lfxsbIy09WpY4eW+7Zt+qdKiot1/gUXJNUcHTNnzpzJr/tv1zOfvKD//myRjJgUjySCQ7Yhw5TsmCGHO9EUTStRDtdEN6mXt7MGtz896f/tTqdT1Y1BvfLc0zrvR2PkdLmU0zFX77yzXtU7tx302Csm/1JDCs846DaP16sn//pHdelzunoW9JLX65UjJVPrX1/RKkyu/9VtSk1Na7PvozVvvSWn06V22dmq37JZ9Y8XKTMjVekZaUpN8cvv98vrcctpGgpvfF+RuK1Yhxy5PR7V19erb9/D7077/HPPquup3ZSWlqadL69U49LFykr1KTMjXampKfL7vHK7XbIsh4IbP1AkO0exlBRFYzH5/X517do1Kef68ssvq+T5vysrI11ul0upKX5Z8Yhee3mFSt5YrfxevZSdfWwvDbN50ybdNv1Wvf7CMqW6HUpLTZHb5ZLX49HOz7YpK6ezevbsyS9P4ASwckeJ1tSWyRFL5Lcdl+yYZDqMxJ9WohSaliTLUMwtDUrrqTNzhrWJf39e165adP8ftXPvPg0YPERp6ekaft4F6lzQXyntcvXppveU3rmnbrn9d/L7D96om5qWJvmy9Jc7fq3804eqc14X9es/QN36FyqtQ552V+5UsKFGP/v1XTpr+Ig2/T5avHCh8gsKZMTjqnr0YaVbhjIz0pSelqKUVL/8Po/cbpdcTqdM01R460cK+1LkzO2kF59/TpddfvmXHp4RiUQ0b+4cnTNihKo3bFBg2ZPKbpeprMx0paWlKsXvk9fjkcvllOkwZdu2Atu2KZzfS06XS598/LHOOvvspJtpLBbTHTNnKMXt0IEHWhmGIZ/Xo3iwQU8uWig5XerVu/c33vMnHA5r0cKFmj3zVqU6Dfn9vlaHeHncLm14v0yXXvbjpDqchkJImBAmhAlhQpgAOISacI2e2bFcGwPbZAYN2fHmTwKVOClcPCo53IYsryFHiiXbthXzSN187fW9dkPlcXiSfgapaWkaMOxsPfrnuXr68SI1RaLas3uXdlVW6uW/L1S4sVa/mvOA+vUfcMjn9zntNNWH4rr31snavHWbGhobFQg0qXzjB9q6YZV+POk3uua6fzvml0w6kcTjcT1w330aPHSomiq2y16zSpmZ6cpIT23Z4JrIWJcsh0OmkcjZYHm5gqd2157aOp0+aJAyv+Say9u2bVNFRYVy0zO07y/z1S4jTVmZGUpPT1Nqik9er0cut1MOh5k4U3ssrlBNjZra5yjm8+n14lc1fsJPk+78FZWVlVq04H5lpCUOSwqFwvp85y4FgyG5nE45nZbS01L0/rvv6N333teQoV9/b6aqqirNuv02rStZoZz27RLrWeGI9lRVq6q6Rl63W5blkGEY2lGxXaMu/JHS0pLnQwp2Gf0OwmRHZKdsy5Yd/+KTQNmJILHjkkzJ4THk8DsUj8RkO6TK6C7VhGuV4UpP+hl06nyK7n5oiW6bfK3WvPyCLr3q58pu30G7Kiv15opnFY9FdePtv1f7Dh0O+fx/ufpnqqup1m+uukTfu3i8zjgrcbza7sodqtn5iX486Tf6l6uubtPvo3g8rtJ//lODhgxR4/ZtclbvVUqHdkpN8SslxS+32yWHaSoajcntCso0zMQJYVa9pqbMLHl9fm3fvl3du3c/5PIrKip0avfuCu7bp9olj6t9RroyMxJl0Ot1y+FwKG7HFQqFZEuKhCNq2l2lpsodMjKz9OaqN2TbdtKFyc6dO7V7+1Z1PaVTS5js3LNXbpdLWRnpcrmc6tK5o55b9LDe27BBv57+W7Vr1+5rh8ldd/xOn215X1075zYXxIj21dQqFA6rY/tsud0uGYah3du3aufOnercuTO/RIHj6NOG7XoitFoyJcMpGXFDdsyWaTUf/pHtkMvjkOk0ZEcSj4mbUT0efE0Tmya0iQyXpP4DB+pvy5br5eUv6tnFj2nHpneUX3Cqxpzu0qJXD/9ct9utG26eqo82bdSWt5Zry5vPq7amXnmn/UBzlqzQoCFD2/zGr4aGBnm8HoXDYYV27JDH45bX65HP65XP65HHncjZWCwmy+GQZCgaiSoQCKqpbKPSMrNUcZgM3759mzp0yFH1m6uV4rKUnpqitNQUpab45fG4ZJqJdQKzeT0h0BSUy+2SWVWlsD9F8XhcDfX1iY3wyfTz+cknSktJkZQ4eeOe+qCmzrpTgaaAXn15pTa9u0an5OYoPS1Vez7dohuu+4X+MHeeTj311K/093zyySf6zdSb5bFDSk9LVSwW02eVu9RnyDBdfs0P5fV59d933alTXE4ZhqG0FL8+/fSTpMpwCiFhQpgQJoQJYQLgEE5N6aKAGZXbYcg2m88QbkgypViOKUc7S2a9ISMu2YYt02XKdhhqdMTU1d+lTc0iu317jb/qaqVlZOjPv75O9/7uFwpvXa0V6yq1e+fOwz43Eolo89sl+s9/u1jn5KfqV/cs1uALLm51mIgk1dfXaeeOSlXv26dQKCjLciotPU05HXPVLjs7Kc/CHovFmj9RCitWVyvLspo/wXLK7XLJ43HLYToSOWuYisXiCvm8SvH7Vf/RFvnO+IGqqqq+dPlVe6vkdZiKrHpd2Vnp8vt88vu+WD+QIcXiMcXjcQWtxCdnTsuSPxbR7nBYMgxFY7Gkm2tFRYU8nsTx9nX1Dbr0/12hLl0Sh6/06t1bWzaP1l/n36c0p+T3eWXFopry84ma/ef5ysvLU3V1terq6hQIBBSJRCQlDnXyer1KS0tTZmamKioq9Kt/v0Eds1LkcLjU2BRQXUS6Ydpt6nXAsYljxl2h4meeVHpaqjwejyq2V0hnJs8sKYSECWFCmBAmhAmAQ9jUUK64Q4qbUtySzLDkcBoKZ0nRrg45qhMnmokHbRmWFHNKpt8h24xrS1O5zkgf2uZmkpmZpVAwpMZASB5fpgb0yNBHWzYd9jmVOz5XJFCn7qd0UDxcpQ8/a9DFp3xxzeVIOKz1767TyuefVfGTDykeCcq2JdOQ4vYXu+n2PftiXTzuJzrznOFHPFHaicTlcikYCikSiciOxRIfCBiGHA5TDodDlumQw2FJpi3FpYgrmtjw63HLbZry7NmtQCDwpcsPBAIKbdsun2x5PB55PR65XW65nYnj/iVDhiSHwyGHaco0zMQZ72NRhcNhRcLhb+0Eat+mHTt2yOVMVJnaugb17Jl/0P29evfRzN/fpYcW/E2VH5UpOytTp3TI1C8nXyfJbj43gluWZclsfpPFbTtxlvVQSI1NAUmG8nLbS5L27qtWbs++uvnnvzjo7PmS1DM/X8vqGpSeliqX09KOHTuSapYUQsKEMCFMCBPCBMAhbAttl+1QIscNWw7LUExSsJtDVrg5W2w7saHXZSjqtBV1J3JnW6hCZ6jtZXhu586ybVs7dlUrPytLBZ1T9fjKFxS+baZcLpfq6mpV+fnnqqqqUjAQkGVZ+vijj2SapnI7ZKj6o62KROPKyekoSSr95wbdN/sula9ZocIBbt38U5+65aYrM80hl2UoFpfqGmP6fE9E724u1rxbXtADWV01afp/adT5FyRF9vj9fsVjMYXDYdlOl2KxuOy4nThsyDZaMlaGIdM0ZTkcshyWnFZiw6+1c6dSmvdmOZTU1FRVfViuji6XXM0bbK3mvDZajpA3Eocrxe3EBeujMVWHExmekpra6hwOyWDPrl0te4g1NDa2OqQjFAqprGyj9u3do2AoIsvpktNpqXd+dxmmIdMwJcNo9UGBbduSbStuJ16neDymSCSqYCiiqj17VFa2UQNPH3TQ2cCz2rVTQ2OjJMk0Te3ds5tCCMKEMCFMCJOTK0yAtqi3ryCR3w5bMachM2grnmIo5pGscOKk4XFLiodsySnFPLZiHkPxqNTLl9jAFIlHZBiG9oT2KsfdQYZhqDZcl7SHhHQ+5RR16DFQa0u3qs+PCpXfKVX1uzfouaf/rrL332vZKGsYhizLSuRFLKaBp/dWmt+jj/bukSSlpafr0Qf/pgfvmKqzhnh0469zdGquSwf+OjXTfyhH3qXKDNer86d/U2Gfj3XFqJiWr6nSPTddrdVjfqabp9+urK95bPd36Yxhw1RbUyNfSorC4Ygikaiikahi0ajisbhs05FYH7QTeWsYhkzDkOVwqPHzz1VeXq6iBx/Ujsod2rsnMcPs9u3VKbeT9uzdq5rSfyq/c6eWDbaGjMRlUuxEJsWiMUWb/85wOKJwJKIGy6n6+nqdM2JEUr4X9+2raslfw2HJ1Zyp8Xhc6999V4sefkg+h6322Vnq0ONUmaZDpmm2fBmGIaN5HWD/cmw7cZ1wOx6XbSfWd+JxU6bpUH6PUxWJRPTswoe15NFH9C8/u0aDhwyRaZqJPDcdLcvad5i9siiEhAlhQpgQJoRJmwwToC36OPyp4g7Jdkq2lfizoZcpxZu3kRmJjb1xr2S7Ext0o05b8Zj0P3v/ovzqbqqPNKi4aZ3a2Wn6kfdcra5fpy6OThrsHyCn4ZQlh0Z3PVdOhzMpZuJ0OjXumuv08N2/0uXnn6Hu3U5V+4wy/fk316uw8DRNv+n/qXuXDspKT5HTcihu22oKhORyWrJjEe2rbZAk3T/vHpW+8pRuvjpTIwanyPF/TgVgWLlyD7tDhjNxVkgr7ywF35qqLL2jCedlaGgfr+5++DH98vptuvPPDyinY8cTem7Dhn1fTzyxRD27dFUoGlcgEFQwGFLI45HTisho3vwai8YUi8Zkx+La1RTQM/uq1WXClerdu7dyOuZo2Pe/L8M0FWveSFy9b5+cLpeqLhmjZSte0uVpqWoXS6w3RaOx5o2TtsKRiEKhiIKhkILBsILhqPZ5vNq1Y4euvjo5T8wXaGqUt/n79HbtZZqm6mpr9UjRQ9q2qVR5uTlyOp1yOBwyTUfiT4fjgAxvznHDaLlW+P51nsRXvDnD44rHYorFEo/tmNNekUhESxbcr7dWD9DPrrlGKampSs/64hJVwaYmCiEIE8KEMCFMTq4wAdqigf7TFHdIMVOKGbYMp62w1yFXU+LyUbYlRaK2zFRTUY+tmEuKOW0ZQenx2ndkRt+RIyqZUcmM7tab0Y9kxiQztkFm9HmZMUOdlKZ4MK4hnQepXUpWUmT5iNE/1IK7vVqxqlSXndlTt/5itKz23dS3Z+eW3ecP5PcmNrbZ4bD6de+k73er1VsrntJNV2Zo1NAv2XvF1VmG0/tFpnsz5Rr0W4VWXypJKshz6+4pHfRfC1bp9ptv1Oz7/3pCHwoy8PTT9T/z5qlT51Pk6tZTjbs+U6PPJ7fLlTgzeGz/OQBiCoXCKtlWoU9PL9QNs+5URkaGYrFYywby/blimqYyMjOVlp6u/IICVY8cpReee1a7Pt2q83v1TCzP4VAsHlMwFFJTU0CNDQE1NDRqX1Z71UejsuMxnT5wYPL/sNq2dlZW6o+z/6A0l6G8Th3lsCw5HFbicJfm/HYcmOEtOd6yn5Vs2Ylt63ZcdvyLDI81Z30sltgrS5LyOnVU7ecf647bb9N//Gqa7MQW+aTEdQi/ZS7DqXtrl8qMSGYk8YZt6uKQFTZkhgyZMhRrsmWkmop6pKhHinhsKRzXB02f6/X6jVrb8KH2BhpV0VitV6vXa2vDTm2o/VAv7lql53e+rnW71itPOUpzpcrtTJxc5ETXIaejlj76iNLTPRrQt6d6nZKqH405X/9y6ZnqntdBGWl+uV1OOS2HXE5Lfq9bHrdTisbki1ZrR0VYa99do//4aYZGF6bKPMR5YgxPbzl7jmv5wTWcXpkZAxSrWCRJapdu6cwBXr2+erNee3Ozzhl9njyeE/eaUdnZ2Xr8kUfVvlMnOUJhuaurWq5ZZBimZNuKxeKJXRODIb3y8afaOGCIrp7+W31v2DBlt28vj8eb+FS6+ReclDiLa7vsbPXu00fdzjxLL1buUmjbpzo1I12mYSgetxWJRBQMhtTY2KTa2gZV19SqMjVTn6WmqrG+TjdMmfKNr9F3PDy99Ck57agkqT4Q1oBBgzX7jv+S3VijDtlZspzO5jCxZFmWHFbzn82h4rAczfebcpgOOUyHTIfZ/Eli4hPaLwrj/l1Lm3fxNQyl+n2q27tLK18uVr+BA/XmG68pzZdYcYqZTl02diy/RIHvWDAekmUkfp95TLfW1JVoW6RWjpAtw5RCuQ45oobMqCEjKtlxQ3HvF/kdky3XPltWk+QMGbIChqwmQ86QEl/BxG1m0JAZkRrCIT1X+Zr+smWxPv9sq87PH33Cn/zM7/fLnd5O/zvvj7rgvLPU099ZHbq2k2Ee5v87Epfx+m6lBLP0fV8Xle7eoy17anTuEL8cjkM8L7JDcnSXI7N7y4kATE+GottfkWKJPSh8HlODe3u17Pn3tbPG1vfPPueEnZ1lWWqXna2S4leVcsopcm3ZLMs05TBNyUjkcjQaUygU0dKyLfJcMUGXjR0rt9utaDTavJH28H9alqVeffpokwxtXbtWPdLTFY1GW8pgXX2DamprVV1br81du+mTHTt00823qGd+flL+rK5Yvlx2OLHxdNu27Xpr1RvqmJkiv8/fnNfNGX1gfluWrP23ORKPMQ8ojKYj8ZoYxhcbfnVAjh+4G5otQ07LktOI6blnnlVDTbU6ZGcl1q3SMnXxJZdQCAkTwoQwIUwIk5MnTIBkE7fj+qjpY0nSR8Gtmlt1nz4IbNIDdQ/r7aa1WhvcoNpYrSqjO/VBaIfMeOJs4ZFsS46Q5AhLCtuST4p7DcXctqI+yVkdk6tWcoYNWcFEbltBQ86wITNsSOHmTxejtuyYFAsnLnJvx6QPGz9TL8cp6pnT44TP8eKVK1S54V1dntJN3k+CkmFLOd4vf4JpSPuCMj4LypXiUo8OKXp4zSZ172KpS86hj+WP712paOX7skMRxRv3KrptheL7XmhenkuyY0rxmup+ilP/+5dX1Xvo2crrcuKeoT2vSxe9X1qqnTU1cvtT5dzxeeKwg3jibOKRSETLt36i1AlXafiIcxWLxb7yVzweV8fcXH1sOrV9/bvq5PEqEAiqobFRNbV12ldTq0875GpzPK4hQ4bop1ddlbSX7lq9epUCtYn1uYbGRuV1bC+vxyOH5ZBlOWRZzuaNt4fL8QPzu/nEec2fHpqGKZmHzu/Exl1btmw5TFM+t1OBYFCZGemyJWXlnqJRo3+YNLNkl9FvGCZbmz5RljNTn0d26LG6J5XryNG6aKk6O3KU4UhXH2e+Ors7yLa2Ke5OXIMwsfYuKW4rHrOlVENxTyJMYh7JuTcuZ5MSpTGixO4mYUOOuKTm0LBtNS8rcYH7/burPb1rtc4rXaELB5x3wofJrp2VSnd45X5jjxS2pNPCUv+sw7xbTSndKWNrSP50r6acNUTXP/OM3tnUpLMGHPqEJpFNv1S0YpkcuefL8GQrXrNZdqj8izCJh5WTZen/t3fe8XFU58J+zsxsk3bVi+Uqy0Vyk4wLxY1mXMDUJHQwYEwNhBtICCEJcC9JCIRwL4EvlxRMgEACCc1AMMXGHWyKsbEtWbJsS7Zl9a7dnZ2Z8/0xq7WFbUrK90niPL/f/ma1fY9m5jnn3XPe99ZLM7jz1//NzFNmc8KMGb22zWbOmsWaNaupbmzEM3oMnooyt0ag7RD1mxi6zjtVe0m9+DJmnXjSYVNMvux28pSpbLQdXl77LqcOHYJlWYQjEVrb2mluaaMqdyBVnZ0cd/zxzJw1q88ew2kZGTS1N7jBgYCPwTkZ+H0+d1Cnaz2mmnxWHLqmx4VxiCwA4tNFHSmRjoOt2QhbYH/mMe5GIqWG3+djSE4Ge20zPlCEtPR0dZJVKP7FSCR1XfUcCNdS0VnJsy0vs8ze3p2bC6mBIwBduss9NAh26mS1+oj5IZIt0SyJrTl4Wh1kGGS6hqVpbkI4TeLf56BHQNcEmi0QEdzSUvGD2zYPTiuzoyD0bq9LuhyLhavv4v6WOi6bcRGG3ju7adVVVfztkXv5yblnkx4zINMLnZabEe9ogV0BTMqC/CCEvAzZX8+MXdm8sb6V6cXJHK3HIjvWYpWvPex2/7S/Ia0o0Y2LKRkJc2cEeOzB+5g09W89Enb1JnRd55b/+C733PUTthoGelYu1B8gZsZICvip6OgkOn8Bc2bMTARqu938VT1ePHEiKxobSdq2mRHJSXSFI7S0trEvmMJGBMOHDOaW734XXdf77PGcnZ1DfWUpuq6RP2RQz2mh3dd7TBs14jN79B45AcQhyeG6l3x0z6QSTtzf3XuoPLjRpIMef7zP5yV/iFs72LEdsrNz+lRbqgGhkomSiZKJkomSiULxtaDD7OSc9VdR6TQmvO2XrquliK/118GRAul1//a3aGid4O0CYqDFBJEUjVjQxttsIyM6TggsnyC5IoZR77j+1iQabpI4NJAOODHZ/cMCdjQeH7a6M0a7j5ESvrvxUboiYW6Ys7hXBnffW7OatIw0jv/WZIgCQY/r7i/6qLqAbPdXRCMtjxPGZPLAc/W0d9qkJB/ukpjlrufyeg5/YbPsWXyTb0FLmY7T/Cpnz0ph2c/WsOnDDzhuWu8t4hoMBvnPe3/KU0/+kZeee47pHj/DGpqwDIOPiidyw+ln/EPOPtL2+Gkn8OLuXcyt2Y8es9jlD7Cio4OLvvlNLrt8IYFAoE8fz3kD89hsWei6t0fCN01zl3G413U8HvfXQE07GOjVDvG40ETPNYROt8NtHFskcjR0LwWRUqJpDrrmlonrfs/upHIxy2JAXp4aECqZKJkomSiZKJl8vWSiUPRmolaUzfs/5bGtT7CvowmfjGdRjgdupQYOEjwCR5NIDzhRt4awKSxI9yLseGyyEzw6WCEdO8kmNsRAdwTBBoG/TYewjdBAM9zpZGiux6UdT1YRn9UjnYNb2+x2uExkd37o0z9x2cyLSEkK9br2XPH3V5l76hQCfi/8g0vv9UCQvAzXIQ2tBx3eFXXYsLWLFR92UVppYttQMNRg9rHJzChOJsnvTm906p8h/NYbYDcBMHSAl2NLfLz999d6tcMBAoEA11x7HXPnzmPlypWsXf4OW9ev4+fxBG2fDeh+0fZo9zmOZOZpc/j9L+9nwqwTOf7kU3jipJMYlp/fL47rQYMGE45G8fm8bmb1hMNdN58wZRLHTZpIVoY746apuYVN28rYtnNXT49/ZgaPI2XC37YQjC8cxTFjC8lIT3P316Zm3v9oE+s2fBB/L9fjTjy7eCQaZfDgwWpAqGSiZKJkomSiZPL1kolC0ZtZtXMNV2/6cXy5RbzCS9yVIr6SQ5cCJwyGV+AYYAuQhjv1O5pqYuUaCDueHdwArc3Byvfha9fwtIO3S6DbHpxUgVVvgiFBl27/wHLDRU68HyAd93OI+Pp33ePeJjTh3ufAoyf9uFf6OxKJULl5I+ccN+efeh1h+AgF3Iyq4aib5GxXjcn9Tzayt8Zi1lQ/C89KQROC8uoojzzTwtLVHfz4qmxy0uPd17i/cf+FzJyYxEPPPM53vn8HwVCo1++Xw/LzuTw/n2MmTeKdt94iPSPjMH+b0Sgi7ocv8njMNBNZxLtvD4ZClMw7netvuIHxEyb0q+N66LBhdHR0kpYSSpSB6s7+/Y0F85k6sZjWtnY+/GQLAIUjCjht5gnkZmex9uPNPWb7HDrLx+3/uP6eNXUSxUWjDnuds+acyoCcbJ5/+VX3PeOvIaWko6OzV69lVQNCJRMlEyUTJRMlE4Xia0U4FubFytcSztZ0d+mFZhwyjTse09EkyAgYHnc5hmOBg0DUWogmi8hQg1iGjrfexlsaJVrixwoK0AUySWK0SaiPIgIC4RFomnufbYKmC3f2wCFBXSd20Nvu+mMQGjxQfBNzJp7aK9tTOg6O5a5Z/6ccLgR6fImIaUmqa2Pc/nAdI4d6uPOqAQzOPlh+Y97xIWYPqeFXS2M88nwjd1+dy5HyoIwY5MWJNbNvbzWFY8b2mX30L88+w5lnnY1lWT3c/OHGjaSlpZGcnEwoNfVgCYTPONyyLKRtY8VilJaVMrqwqMdjisaM4Zk//Ymf3Xdfvzq2W1pa8Hq9if3J9bjgmAnjmDKxmK1l5Vx76x0cqDkAwIC8Afz2wZ9TXDSKmoZG9tY2oOs64UiU1na3pFlqKEjA78O2BcMGDqC4aBRby8q5/rY72b9/PwC5ubn89qH7mFpSTPnOXXy4eXOPpSNer5eWlpY+1ZYain9YJt1/d68JEiLeoBEwomDEQI+4CWGMWgtPaQQrahFJktBq4d0eRUQdrKAklgbRLImZ7GCFo4gAiCSBHhAYSQLdL1vZEhIAAB4XSURBVPAka3iDAsMv8CQLPEGBZoDuFWgege51BaYZ8KspX1+ZpKdo/ObOAXzvkmxOPyGFeceHuOlbWfziQoNIkyuTeMWFo8gkwr691X1qH/3Ls89wwrRpWJbV4/LBhg1U7dlDQ10dsVjssPu7L7FYDDsWw4rF2LL5k8Rt3dtumfQ3voxMzrj4Si6++kYuvvpGzrj4SraVlVNcNIqCIYMS5SiiZozahiZqG5qImrF4aQqDgiGDEjI585JFidc5/cIr2FpWztSSYiaOH+dOOe3jMlEoeqVvpOQPG//I0ub3ErNnpAO6rzuA4z7OibkBVyndZRiOBYQlxAQyCtISEAMnScNbZ+OptcHQ8FbGsDstYu0RrOouYnu7sD3gBAW2DyyfW39YiyfRFJpA9wq3v6CD4XcHgLoPNK/r9x+NXcil0y/stYnh/IEA2flF7K9rTtxW19jGky+uorW9Z/3UmvoWfv3kMlrauo7QF7AxLSfRGX3kr40MGWDwwyuyewwGEwHQ4UEumBDhw09NahpjR/xsWakGUkJdXV2f2Uf379+PYRj4/P4e7m1qauLhB+4jIzOTcCRMc3PzYW7u3ra2tmI7NplZmfzlj0tobW3p4XjD46Fqz+7EgKa/sKOsjFAwuafDheC4SROJRqM88NgT/Pyhh7n+P24F4EDNAb53z8+IRKOMHVmAbhiYMYvUrAEsvvFmFt94M6nZAzBjFrphMHZkAZFolMeefZH7/vthLl20GIDa2lq+d9dPiZpRjp88scd7A4SCyewoK1MDQiUTJRMlEyUTJZOvl0wUit7I6p1r+UXlk+4KXwl2VGJFJbFOiRWRxDrcrWPibqPx7N2mdI9JQEvSwKcRmejFCWpoHg3h1d1LQEe3BbEhXmwPWEFBLE0QCwliqQIrFC8nlSRxfNJdImG4fQjNI9AM19ua7gZ3/V6D+WNP69X1hIUQzPvG+bz65vt0RbozI0veXrWJ7/7XE+zeWw/A1vK9fPvOx9hVdQDDOLy7aYfbqarvdM+tTRaflsW46qw0goEjd001QyfV7zq/M3LkqK7PKwgmC9r6UEBt27ZtFBaNOSxQW1dbixmJcMMVl3HX7d9DCHHUoC7AT75/GzdccTmWadLY0HjYYwYOHsy2bdv61fG9s6ICn9fTo9avQJCRnkb1/hp27d7DW3//O9s/3ZJ4TllZOftqakkLBdE0jbaOLuYtWEAoJYVQSgrzzlhAW0cXmqaRFgqyr6aWHRUVLHv9NXaWlx/sP5RXsHd/DRnpaW6egEM+g8/rYWdFhRoQKpkomSiZKJkomXy9ZKJQ9DYsx+KeTb/CsV1v26aMexN3lk2ywAiAJ0lgJLsuNQICw+e6VXgAj+tcmSyIpmvurJ+oQNM1hE93L4YGhsBO1zCzNcx0DTNdYIbADEiiAYgmg5kCUb8kFpTEkiV2UOIkSWQSFKcP5fGp/8mH5y1lVO7IXt+2p86dT1NLlOdffw8pJbmZqfzP3YsYPDCb677/CH94bgXfufM3TD9+PPfeeiHBpM8mC5CYDXtY/kktxWO8dIQd0tM0CocePcO31WVS1WIgBGSnGUftX/h9IuG1vkDZ9u1kZWUdFqgNBoN4fW67nXHetzAM44gBXcuy0HWdE+fMA8DweAkEAkd4vRBlpaX96hiv2rM7UUOxO7EbAmzbwe/1sW/Pbp754xLefevNg3uelPh9XmzHwTRjjCuZSCgUSgRlQ6EQ40omYpoxbMfB7/NSvm0rf3piCetXvfuZPqMP23Y4+NbxmteaRtXuXWpAqGSiZKJkomSiZPL1kolC0dswNIOrCi5wA6c+d3aNZrgzbRzrkPV78dwAQuNgwjYNHI97sXSJlSpAgt4lMUyJEQHNxK0VHN8amoaOwEyFaAjMYHwbkkS7L+mSSIokkgqRVJBpOsdkFXDtmMuZUzSb9KS0PtG2AwcN4nsPPsYTT7/OUy+tJhwxSU9N5s4bz+WS82fz9J+Xce2VZ/GdhfNJChzuZbNhL5u3V/De1kbOPTFEbaPF8EEGhn7kmU3Sdthd1sTzW72ccWIS6aEjB72llIQjEo/X22f205qa/QnnWpaFE7PwmzbZSSG+95O7uen2HzJj1qyjBnS7L5OnHsuFV17NpVdfQ14wlXR0sJ0eM30O1NT0q2O8sf6zs7ncoO7e/TUMysvlG+ecedhzzj/vbPJyc6lvaqa5tY3iiccc9pjiicfQ3NpGfVMzebk5fOu8sw97zDfOOZNBebns3V8Td7f4zGer71vnS6WML5bJD7b/TyKBoNDAseOlIOQhMnEOkYl9iEyMo8vEMUF+RiaO44pCSNBsEDZoljz4t3RvE1IgJPgdnUlGPlcPvIg5I2f3mbbtlsnPbrgIIeBb849PyOQv+QN4/MnXuGHxuZw3Zyq6rn2uTO66LpMPtoe/1jIpmTgxkUwG28FvOQmZtLW1MSw//wszj06eeizZ2Tn4/X7ygqloUtJgx7Bs93n9VSbpPnFEmRSPLeIb55zJ315aekSZ7KzeS3NrG/OPIpNPN31EfVMzI4YM4lvnnc1fX3zliDLZvK20X8hEoeiNXHzM+aT6UvjO5p8Rk3YiO7emgW2BrifSACCEu+xDGuB4wfaC7ZPu7J1k0KNxX3eCaJN4NLAjYEUlIibQpECLgLcTIqmSy4Mn8HFsB/ucNiYyhOVWJTd4TuHk5OnsNw9QHBzHpPSJBIy+Wb5n9tx5yEef4Ze3Xcdrb23gigtPY/6sEi49ewbzT5xIVnroSJLFbKhif/nHPPxyGceV+JhUGOCjsgiOPPL7dIYd/vRqA2+u95A3yOCiOUcfNIejks4uSWpaWp9pR03TMWMxvEKjqF2SFQVdaoDG1FA+lWlRtjvhL8wY7tgOZ48qYXwgg4Dmdu8t6VDZ2cxrdZWEY7FEALS/4DgO0N2fc9d1SSQr1r7HmNEj+eF3biAtNYUnn3kOhGDhxedz/cJLsGyLjVu2Mzi/gOycHKSUPRLDZefkMDi/gI1bShk2MI87br6etNQU/vjMc4Bk4UXnc/0VlxKzLFasXe9WBDg0qQjgSEcNCJVMlEyUTJRMlEy+XjJRKHojQggWjJvPkNRB3Pz+j9htNSAkOLabgK07N4BmuOWjbA/YusQJuNctAxw/WJ74L4IRiRGVEHVLQukeiW4JNEuiS4GIgROVeMOCczPn8b+Zv6Q51kJQTyYmLaJOlGxvVr9p2znzT6e25ic8fu9tWDE7cfth/pYSq72RaN1OdldV8eDfShFGjBu+mYOhC/KyDFZ9GMay5WGBXdOS1HXABQtSmHtc8Ig1h7upb7EQAnJzB/SZdhwyZAhNjY3M9GaSG+353X1SMMby047JTjv6uQ4v8oSYmpzTs5MvNEYHM8Fx+O329yksLOwX+157Wxser5f0jEwItyQM3q3R2oYGXnhtGefMP407br6emxYtRAhITkoiYpq8ueZ9Kvce4Npvn5/wdjgcBsDv9yOEYPacuTz2yMO8ueY9Zk87lh/cdB03LbocJCQnJxGORHnhtWXUNjQm3vfQbmh6RibRSATTNAmlpKgBoZKJkomSiZKJkkn/l4lC0ZspGVzMzU1Xcfv2B4kJm3jMC0fEg7gesL0Suzug6wHHF3e6T+L4wB+W6FHQImDF3GRuRECz3Bk8/lrc55sCOyrZ1PEpZ2TNSzjbDd2G+lW71tYe4KkH/5OrFy5gwSnHHPVxZkMVZZ9sYNXWOl5YvZdR+QY/uiyHrFT3H1Ew0Etrm8Peuhj5eT1n6KSHdH50Zc6X+jzl1VF8oSwG9aE6ruMnTGD9+nXowzP47EyRhIvjUz8/L7DrNY6eSNArdBrq6/nmN7/ZL/a71WtWs2tnJXV1dRzYvROvz0tyUoBAIInkpABJSUls+PgT9tXWMW3KJPJyskEIarbvYN1HmzHRuXzRYlLT0igtLeU3jz7K7l2VAOQPL+D6G2+kqKiIyxct5uklj1NaUcn0SSUMyMkCKanZXsa6Dz5i3/4aOjs76erqoiscdi+dXURNkwHDRvD73/2O4SMKOP30M9SAUMlEyUTJRMlEyaT/y0Sh6O2cN+Fs5hbOpqqpmt9sX8LfutZge2RirWAPf3shzWNwim8ct6QvJl1PY0b5VYioBZabTVxooGsCYbolpqTtlpnShcTwwBNtr/Pt2DWkevpnQGfPrl289PxfiHW1csbJk46a2VxaJpHaMu5+egtNbSaXnhninFkpBHwHZ5yMGuIjLVVj7Zauwxz+ZXEkLP+gk9Mu+A6BpKQ+044Tiov5P48+QsmIQqaIAKnyYNBaAlV2J9vMdizH/lyHb2g9QI7hZ2RyBtoh/4u6cDsvV22jvra239QS7nZiRmYmbzz/FB7DgyMdHAm2I2lubWXHzt3UNTZTua+WYCiVtPR0hgwbxknzz6Jo7Fg8Hg9SSh7//e+JxWJMPe54APbv28fjf/g9v7j/AXIHDODm226jdNt2tm/byoqPttLS3ExHeyutTY3s27ef7Mx0/D4voWCQ1FAQXRPEYhbzL7iMCy+6qM+0qRoQKpkomSiZKJkomSgU/R4hBEFfkLF5Y/i292oW2ZdiGAbV0b3UxxqYkjaJjyObSTGCjE0qwqt5GewdhEAQdUx8MQ3LlDgxCVIgbdflmiMw2iVaGIThlpmyTQhHY+yN7CXVM7ZftWNzUxNP/v63vPS/95Gbl821i84hJeg/ulvNCJqAOy8ax51LNlPbaOH5TADS7xVcMDfEH15oZfaUILkZX717Wr4n7GYbv3d+n2rPYDDIBRddxIbNW2gZPZoUR5AiNSzH5oAdpd2OfeGSj+7tc/tLSdU8DPUF0RHs62pjV2cL+6qrueGmmwgGg/1qXywaM4bnOzrJzszA0Aw8hkFXxCQpPYdb7riL0UVFhEKhw84DhzJ+wnjeXfEukUgEgEgkwvEnnJC43+PxMqGkhPHFxT2e197ezo6yUp5a8jix9maS/F5iloV0HJqaWyksKupTbakGhEomSiZKJkomSiYKxdeKwszRB49hDjp2bMqRjztD6FyTfAa/ib2A47juRgg0CcKWiGZ36ji2QDfB8IBlOuR5B/Srdqvas4c7bliMbNvDT3+8iCnjh+P1fL5v9aQUkkbOpNC7iR9caHHPk1vIz2nh7JPTepxPT50c5JWVHSx5tZnbLsk6apK4I9F5oIOnn2ti7MwFTCgp6XPtOnfuPNauWcOB2lpa/H6i0SiWZSWWIHwZf3dvG2SMFtvEtm0sy6K9vY1BgwYyZ+7c/nccFxYSjieFl45DR2cXXdLgPxZfQ3Z29hH7LIeu+Qc459zzKCsto6J8ByAYMWIE55x77mHPOVLfK3dAHgsXLeahB36BbXfhix8LYZs+t8RGv/vuu+9WavjqZCVlMiA5l5xANqNCIyhJnUC2L4vi5HEUBkaRYaSTqqck0sgLINzeyAc125ERt54hIp4C2w9OGPSoBAQaIHW3IP0PhlzTZ5PGHE0mt169kL1b3uWOWy7kxsvnUTJm2FF/HQTQPD6M1IGkamEKcg2WvLqPYMChMN/f43mDsz2s+qSLqlqT48cnoWlfTSb/55lG0ksWsPCa6/pc8pThwwt4442/o+kGnY5FdUcL+zpaiSETouieMvpF2y4rRr0dZV+kg5ZYhNbWFrIyM1l8zbXout6vjuPU1FT+9PTTpAYDICXhSJROR+P6m75Dbm4umVlZR9w3D72toGAEn2zaxM6Kcpqamhg6dChXL16Mz+f7XJl4vV4cKRlVWMTK1WuQloke32cb28PccPMteDwedbJVKHoBmtAoj1awsuoDnDBYUTezuO4RSAu09nhOAQfXPTqYhmSyfwSFKaP6RRs0Nzdz66LLGRjs4r7bL2V0/oAjZgE/EsLw4EnLI8cfwyMiLHmtjqJQmJxsH5rH9YrXIxg5xMMfXmxDGJLxBX7EF2jcjli0ldexcm0TS7cb3PnQYwwcNKhPtWssFqOyshLbcajZv5/Ozk6Sk5MxjHhyt7ibv+jiOA6O4xCLxYjF3ECwaZpkZ2dRWDSGQCBAampqv/D4u++uYO2ataxfu5bde6qoKCujobGJnburueSKRTiOTXIwSGZm1mHOLi0txbKsRLDX5/MxY+ZMwA0SX3f99SQdMkuspqaGfXv3kpWVdZjTD9TUUF1VxcBBg/nbc8/R2t7BvppaUrMH0NbWRnl5OU2NjeQPH64GhAolEyUTJRMlk/4vE4Wiv+MRBusjG6iv78COuA7XPAIZAycm0XS3vBRCoAmBY8DpmdMZk9Y/knH9ackfKF33Gvf/8DIy0r76bBGhaRgpuQxLjlBaVc+azZ1M1GuwWyOYrRHCte0EmhvJ9Ng8udJGYjMyU6Drbo1XgdtPckwbsyVMx54m2nZUs3ufya/XB/jmLfdwxlnnfG6AuTdhmiabNm1i7dq1ZGdnM336dCZOnMjrS5fi8/lwHIeuri5M08RxHKSUnxvQFULg8Xjwer0IIWhva8OMRLj9Bz9g4sSJNDY2snLlSqLRKFlZWX3a5QF/gMKiIo4/4QTmnX46b/79dUYMG8TICZO49fvfZ9So0Ql/J/a/+H7R2NDA1k8/ZdTo0YnbdV1n44b3sWIWk6dM6eHple+uICMzM+HwQwO8mVlZjBw1ignFxZRX7CQz2YstDJY88yzTpk9neEEB6RkZfWJ2lRoQKpkomSiZKJkomSgUii8g15fDtvZtfGjuhCbcQvcekN1F7nXh1iPuPl8YcH7uXEamFfT5797W1spd117KjVfOo7ho6D/8OkLT8CalMiyplefX1uDXNQqCHVgdXdhdYaRlMSTNIjdZ8NRqh4rSDgZ0HYAD9XTubaRzTz2d1fVE6lqwOrvY12bw0Nokxp52MbfcfkefqSFcWVnJO++8Q0FBATNmzCAnx018t3v3bnSPBwFs37qVrOxsAGzbRkqJbdu0NjRQW1lJS10dZiSCJxDAtm08Ho97f2sr+6qrmDR5MqNGjyY5OZns7Gxyc3MZN24c4XCYd955h0AgQEZGRp/cH4PBIIFAAMMwSE5OZuDQYSx98QVmnDqHoqIiSktLWf7O2xQWFvXoqwghyMjI4I3XX8MwDLJzcojFYqxauZKO9nYQ0NTYRN7AgUgp+fijj/hk08fMP/0MhBA9/G2aJn/587NIKQmFQtTW1rJh7Rpuv+enjBs/HsMw8Pv9fcbfakCoZKJkomSiZKJkolAovgTrmt5nQ7QMv0fgMQTSAcdyL0IIHFuiGQItHoicHBrNpJyJff577yyvYOlTv+Hay+eTGvrnkq5pXj/p6ZnIrlr+/F4bOjpBn6QtotEa1WiJaAS9khQvrNxjsGqXF6/mFnpu6dJoi2h0xjR2Nnp4eF0ShSdfwI9+9guCob6Rhf2dd94hGo1y2mmnkZWVhZSSsrIyli9fTl5eHrNmzWLK1KnU1dexb+8+WpqbCSQlJWbz7N22jd1Ll9JSWkrYNEnPz8e2bcLhMO1tbYSCQSZPmcpVixYxbtw4WlpaWL58Obquk5mZSWZmJkVFRezYsYMdO3ZQUND3+5jJycls+PBjNqxfR11dHekZ6Uw99jhSUlIOC/JrmsbkKVMo3baN5/78LO+9t5784cM5fcECiotL2L9/H3959lnWrFpFRmYm559/AYZhHLb0Q9M0cnMHsLe6mldefpk3Xl3KoIKRXHLppST3QW8LeaTFLYp/G3eX3cev972MrwF0U2DHJHbUXVOoewUSieEX7nU/3FN4LYvGL+zz33vbp59y41kzePLR7zEkL/Offr1YSy1P/fllnnqzknPHxJg0yHR//cPNFiqlYNN+D0vLPHh1OG+syYhMCwFoQmLoUNeh84cP/Iw99Xx+/PP7+0wx+u7B2LHHHothGDiOw44dO9i0aRMlJSWMGjUKTdN44okl7KzYSSQSISs7OzEldNdHH1H9xhvuSbSoiNGnnZZYZxiNREhLS2PkqFEsXLgQKSUVFRWJ1y4sLETTNCzLYsOGDYTDYU499dQ+v3/W1dVx+223Ub1nN7NOOpmpxx1L0Zix5ObmHiYTIQSxWIx3ly9nw/vvoek6M2bOYtr06QCsX7eW1atW4dg2U487npNPPhmP13uYTKSU1NXVUbp9Ox9s3Mi7b7/F4GHDuP+XD5KTm6tOlgpFbzxXhOuZ/v75dDVF0LoEtilxYmBFpFtKCjewq/tA9wvuHns9V024rO8PhFev4kdXnM1LT/yYlNC/Jgt3U90Brvz+b2lr6/qnXueCW+7m/EsuIyMzs0+0ZVtbGynx2rLNzc0sX76cESNGMD4eCIxEIrz77rukpqZSWFjIqpUrefrJJ0kOBsnIzMSj6xwoLycWiZCePwxLQmNDPY5lccWiq5k5axZlZWW0trZy0kkn4ff7sW2bLVu2sHPnTk455RTS09MP+yx9CSklO3fuZMuWzWzZvJkRI0ZSUlLizpjKymLO3Llo8WDu0WZ9fdnZYEcbJnXfbts2b735JnV1dcycOZPNmzdTVlZKcXEJE4qLGTlyZJ+YeaYGhEomSiZKJkomSiYKheJL8l+b7uc3u19C7wQrLJGO63ChxZeCGO7sn2FJ6az+1isYWt9P6L5+zRruXHgmLyz5MWkp/7qyTDV1zbS0hxECBMLdCvfXVc298eD5T0qklDhSYtsOEdNix64anvnrcmxPOrc/8DDTZszsM+fLHTt28OmnnzJ79uyER+vq6lixYgUnnngiAwYczFBrmibl5eVs2byZT7dsprq6GiQMHjqE8eMnMKG4mNGjR+M9ZJbTgQMHWLlyJSeffHJiFlFbWxtvv/0248ePZ/To0X16n1y9ejUFBQXk5eUlEgFKKVm/bh3Llr3BRRddzNBhw76Urz97/xcNjbrvr9qzhz8/+yxz5s7lhGnTEq/jOA4HDhxgV2Ul02fM6BPtqQaESiZKJkomSiZKJgqF4kuyonol3/3kXhpaO7GjEtt08wDYpkTzCHSPQPfBGVlT+d1ZD/eL71y1Zw8LTy7hnjuuYNbUooRTwc2YDv//vNna2s7jv1vCK+v2c809/8P5F1/a69e4v//++5imybRp0xKftbKyks2bNzN37lwCgcCXcsgXuSkcDrNs2TKKi4sTU0Nt22bdunV4vV6OO+64fnmMtra28urSV6iuquK0OXMZOWpUj/7NP4NpmlSUl/PWm8sYPGQoC848k7Q+MsNMDQiVTJRMlEyUTJRMFArFv4CIFaW2vZbz3ryamq521+GmjM/0EWg66D7BXWMWcd30q/vFd5ZS8uhDD/LS//6cK86bxtShBgNCAk2TbrEsTQdNR2gGQvPEr+ugGQjNSPx98LqWeDyahhBa/PHx+4SW2B7hwyClDbaNdGzscBvhqg95a9MBHnmpnCvvfIBLr1zUq8tHVVVVMXTowXwK5eXlVFRUMGfOnH95/8O2bd58801GjBjRI5D72c/QH2lra+WDjR+wauVKsrOzGDtuHIOHDCUrKwufz/eFfSApJdFolIaGBvZWV7F161YaGxqZOWsWk6dMITU1td+0lRoQKpkomSiZKJkomSgUiq/Ii5tf4dsf/xw7ArYlsbrcoK7Q4e5xV/LtU67tV9/XNE1ee/klljz4M9rrdpOS7OGYkWkMykoiPegl2Wfg9+oEfDq6JtA1gaHFs3wLvlJt4EO6qYmgsZuYyzkYTD7shAuvf7Cf197bz62/WsKCc87tE+26e/dutm/fzty5c/9t/Q7HcVi2bBljxowhPz//a3m8tra2smfPHnZVVrJn925aWppxbBt/IAm/34dhuHV/LStGJBIlEu5CNwxSU9MYlp/P8IIChg0b1m+9rQaESiZKJkomSiZKJgqF4ivSEelgR20Fz3z8V/5Y/TZ2FBwHhgdSeeeKZ8kIZfTL793Z2UnFjjIqdpRRWV7OrrLtVG7eSLS9odd8Rs3j54UNZb3+fNvZ2cmyZcs466yzErWD/11YlsUrr7zC3LlzSU5OVgdwvG/jOA6maWLbNgC6ruP1etE0rVf/MKAGhEomSiZKJkomSiYKhaKXYNkWr3+8jFtXP0B9LMyfZ9/LvEmnfW2+v5QSx3GIRqNEIxFM08SyrHidWwsZX78vnf933U2fz8fQPhC87C4H9e/296Ee13VdJSpTqAGhkomSiZKJkolCoVD8q6lpOkDp3jJmjZ+BrumqQRQKhRoQKpRMFAqFQqFQKBQKhRoQKhQKhUKhUCgUCoWiF6IWuCgUCoVCoVAoFAqFGhAqFAqFQqFQKBQKhUINCBUKhUKhUCgUCoVCoQaECoVCoVAoFAqFQqFQA0KFQqFQKBQKhUKhUPQj/i/X01YoCV0rOAAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![blue_red.png](attachment:blue_red.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We expect the blue ecosystem to be more diverse than the red one, since the proportion of species is more balanced there. We also expect that the diversity increases when there are more species. This is exactly how $S$ work. The more species there are and the more balanced their proportion is, the greater $S$ is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### PE and barcodes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PE is just the Shannon entropy of the barcode. Consider the barcode $\\{[a_i, b_i)\\}_{i=1\\ldots,n}$ and calculate its bars' length $l_i = b_i - a_i$ and the total length $L = \\sum_{i=1\\ldots,n} l_i $. Then,\n",
    "denote $p_i = l_i/L$ is a finite probability distribution and"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ PE = \\sum_{i=1\\ldots,n} -p_i log(p_i) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the same way that the Shannon index, this value will be greater if there are many bars and their length is balanced. Note that it does not keep information about the birth or death of the bars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PE in Gudhi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are using PE to finding differences between two family of point clouds: one following a normal distribution and the other following a uniform distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAipUlEQVR4nO3df5CV1X0/8Pdnf9xddVH5wq4QYFm0iOyyCuyWoB0n2JiEtib+mKL2S0IoaRyqTUqqrVUTCc3UjlOniTSSHSb+iBmaAtqOnbQ135iRr+mgDQuu8ksDdFhZxO5KQUFgd+/u6R9373rv5f54fp9znuf9msmQ3b0+99zn3nue83zO53yOKKVARET2qtLdACIi8ocdORGR5diRExFZjh05EZHl2JETEVmuRseTTp48WbW0tOh4aiIia+3cufN9pVRj4e+1dOQtLS3o7u7W8dRERNYSkd5iv2dohYjIcuzIiYgsx46ciMhyWmLkxQwPD6Ovrw/nzp3T3RSj1dfXY/r06aitrdXdFCIyhDEdeV9fHyZMmICWlhaIiO7mGEkphePHj6Ovrw+zZs3S3RwiMoQxoZVz585h0qRJ7MTLEBFMmjSJdy1ElMeYETkAduIO8BxRJQOnBnH3pp3Yd+xDtE69GBuWd6BxQp3uZlGIjBmRE1Ew7t60E7veOYGPBkew650TuHvTTt1NopCxIzdIS0sL3n//fd3NIMvtO/YhRkYz/39kNPMzxRs78oCk02ndTSACALROvRjVY9/s6qrMzxRv7MhzHD58GHPnzsVXv/pVtLW14bOf/SzOnj2Lnp4eLF68GFdffTVuvfVWnDhxAgCwZMkSrFmzBp2dnXj88cexZMkSfOMb30BnZyfmzp2LHTt24LbbbsPs2bPxzW9+c/x5brnlFnR0dKCtrQ0bN27U9XITZ+DUIJZ1bUfb2hexrGs7Bk4N6m5SKDYs78DC5om4qK4aC5snYsPyDt1NopBZ25GH9aU8cOAA7rnnHuzduxeXXnopnn/+eaxYsQKPPvoo3nzzTbS3t2PdunXjjx8aGkJ3dzfuvfdeAEAqlUJ3dzdWr16Nm2++GU888QT27NmDZ555BsePHwcAPPXUU9i5cye6u7uxfv368d9TuJISO26cUIetq6/D3nVLsXX1dZzoTABrO/KwvpSzZs3C/PnzAQAdHR04dOgQTp48iU996lMAgC9/+ct45ZVXxh9/xx135P33X/jCFwAA7e3taGtrw9SpU1FXV4fLL78cR44cAQCsX78e11xzDRYvXowjR47gwIEDgbSdymPsmOLKd0cuIjNE5GUR2Scie0XkT4NoWCVhfSnr6j4evVRXV+PkyZNlH3/RRRcV/e+rqqryjlVVVYV0Oo1t27bhpZdewquvvoo33ngDCxYsYF54RBg7prgKYkSeBnCvUqoVwGIA94hIawDHLSuqL+Ull1yCiRMn4pe//CUA4Mc//vH46NyLDz74ABMnTsSFF16It956C6+99lpQTaUKGDumuPK9IEgpdQzAsbH/f0pE9gOYBmCf32OXs2F5x3mLHsLyox/9CKtXr8aZM2dw+eWX4+mnn/Z8rKVLl6Krqwtz587FnDlzsHjx4gBbSuVkY8dEcSNKqeAOJtIC4BUA85RSHxb87S4AdwFAc3NzR29vfn30/fv3Y+7cuYG1Jc54rsiT0/3AlhXAe7uBKe3A7c8CDU3nPSz0laEO20HnE5GdSqnOwt8HNtkpIg0AngewprATBwCl1EalVKdSqrOx8bydiogobFtWAEd2AEOnM/9uWVH0YaFn9zhsBzkXSEcuIrXIdOKblFL/FMQxiShg7+0G1NjCNZXO/FxE6Nk9DttBzgWRtSIAngSwXyn1d/6bREShmNIOyNi0mNRkfi4i9EQCh+0g54IYkf8WgC8B+G0R6Rn73+8GcFwiCtLtzwIzfhNINWT+vf3Zog8LPbvHYTvIuSCyVv4DAGurEpmuoQlY9WLFh4We3eOwHeSctSs7iYgogx25D93d3fj6178OABgcHMSNN96I+fPnY/PmzZpbZrHT/cBTS4FHpmX+Pd2vu0XkBN83rYzaIcg2nZ2d6OzMpHS+/vrrAICenh7H//3IyAiqq6vDaJq9sqlpKv1xahpvw83H900rjshzHD58GPPmzRv/+bHHHsO3v/1tLFmyBPfffz8WLVqEK6+8cny5/rZt23DTTTehv78fX/ziF7Fjxw7Mnz8fhw4dwi9+8QssWLAA7e3tWLVqFQYHM9UZW1pacP/992PhwoXYunUrWlpa8MADD2D+/Pno7OzErl278LnPfQ5XXHEFurq6tJwHrZiaZie+b1rZ25FHfCuXTqfxq1/9Ct/73vfyytgCQFNTE374wx/i+uuvR09PD6ZNm4aVK1di8+bN2L17N9LpNH7wgx+MP37SpEnYtWsX7rzzTgBAc3Mzenp6cP3112PlypV47rnn8Nprr2Ht2rWhviYjMTXNTnzftLK3I494ddhtt90GIFPa9vDhw2Uf+/bbb2PWrFm48sorAbgrffvJT34SEyZMQGNjI+rq6ipWX3TDio0VmJpmJ75vWtkbIw/hVq6mpgajo6PjP+eWl82WpK2urva9rZvb0rdByS69HhnF+NJr44pIMTXNTnzftLJ3RB7Crdxll12G/v5+HD9+HIODg/jpT3/q6Thz5szB4cOHcfDgQQD+S98GhRsrUB5mmsSGvR15CLdytbW1ePjhh7Fo0SJ85jOfwVVXXeXpOPX19Xj66aexbNkytLe3o6qqCqtXr/bdPr+4sQLlMaF4lZOLCS84FQVaxtapzs5O1d3dnfc7lmZ1zuu5Cr08KdnlkWmZTjwr1QA8eDSYYzstVfvU0o/TFqUmMygrDNE4eUxClCpja2+MnFzjxgqUZ0p7fgcZZKaJ07xyJ3NdTG2syN7QChH5E2amidPO18lcF1MbKzKqI9cR5rENzxEFJptp8uDRzL9B7tLjtPN1cjFhamNFxoRW6uvrcfz4cUyaNAmZEudUSCmF48ePo76+XndTyEKRzpHc/uz5MfJinKQtMrWxImMmO4eHh9HX15eXu03nq6+vx/Tp01FbW+vvQNw3MXGWdW0fX0dQXQUsbJ6YP2fCz4TxjJ/srK2txaxZs3Q3IzlY5ChxKq4j4GfCWkbFyClCzARInIrrCPiZsBY7clsEvSiCmQCJU3ELt8LPROMcLsSxhDExcirjdD/w9wuBwVNjv6gGmhf5u+1lPJQKFX4mRoaBd3u4EMcgxsfIqYwtK3I6cQAY8X/by0yAypJ2sSv8TDwyjaEWSzC0YoNiXyCGQsJnQi0SnRh+swY7chtMaQckZ0u4uglcFBGFpE/+cSGONRhasUGxxRVxvsU3hYNaJLEuRMbwmzU42UlUioMYecVFNkQB4mSnWzGb6Ir1yDEsDkak3KyDTGB1jDzUPShjNtGV3ebto8GR8W3e/LJiD9CQcbMOMoHVHXkYndO4mE10hTFyDPX8+xHhjjIVF9lQcLhTUElWh1ZCva0Ns+h+VoThm9apF+fFcouOHF22x9iwQoQ1Q7hZR4RYC6Ykq0fkod7WRpF6FWH4xtHI0WV7jA0rxOxuisbwfS3J6hH5huUd503gBSaK1KuCD+boe7txR9f2UCYkHY0cXX5RQj3/fkxphzqyA6LSGFZVOKBmovHUICd3bRfFXbKlmH6oU8Gmsvtq5uDzHz2kL5UtLpvcnu7HvsdvQfPQIewdnYmvjazBzOYWhkBsVyT0N6AuSVQ2FtMPTVSw0Gd170q9MWenu7qYKueLfmZoOpYMfhfv4xIAwEemxO/JuyJ3yXfn5PFnJ9yTeMFmR65TwQfzsq7tOFppQjLC9lgnZzJsvhzAhtT3cPvQWrPi9xQoYyfcI2b1ZGfcMJXNp5wYfw1GMK/6HSvPJfPznTN2wj1igcTIReQpADcB6FdKzav0eN8x8pituqSAxCTGz2X/ziVtxXLYMfJnAHwfQDRBVd35pEm/kJj6+m2P8Y9huMA55vFnBNKRK6VeEZGWII7liO580sILyU/uBKpT5nVsYTH19dse4x/jaPEWUY7IYuQicpeIdItI98DAgL+D6S54X3gheff1WNVlqSjprz9knCshtyLryJVSG5VSnUqpzsbGRn8H013wvvBCAknWirM4vn6D6nhkwwX//4/bsO74fbjgsWbs++vr8P57R9wfzKDXReGxM2slewv94NHMv1HfxhdeSD6xwPotsVxlSsTw9ZtY7bL/yTswe2g/GuQcZg/tR/+Td7g/iIGvi4LHPHIvCmOxxSb/LJOtZOhoYYVlr79SZsPAqUFcdKQHFxp2V9E8dAi1kpn1rJVRNA8dcn8Q3fNJFIlAOnIR+QmAJQAmi0gfgLVKqSeDOLYVXEyymZou5StTwvBJxkoXqbs37cSfjzRjgRxArYwijWrUGHBX8U7qCswe2o9aGcWwqsI7qSvQ6vYgrE+SCIGEVpRSf6CUmqqUqlVKTY+kE7c09mdqDe8wF1boXuBS6SK179iHuHtoDXaNzsZpVY8eNduIu4qmr2zGgdRcnFb1OJCai6avbHZ/EN3zSRQJe0MrunPJPTI1RzjMSoauwjYhqJTOl/n7CO4YXvvxAhwD0kcnT5mByQ9tBwD3I/Esw++WKBj2duSWxv5MzREOc2GF7otXpYuUseV4o2LqAi9yzN6O3NLYXxI7Dd0Xr2IXKVPnKrTwe3fLC4F29tYj54fHGiZ2mjbVMwn9/D0yLZOemJVqyKT2OuWxxo2JnwvTxa8eOWN/1jCxHobucI8boc8x+L279Rjm9P26OJgbZ+eCIAqF7uySKNlU/jT0i47fzBaPJTN8vy4udhrHjpzGmZoaGYbv3DwPF9RmOp8LamvwnZsrVl/WJvSLjt+V0h4vBL5fl6UJD2FgR07jbAo3+PWtF/bg7HCmEzg7nMa3XtijuUWlGV9Ey+OFwPfr0l08zyD2xsgpcLqzS/zKnTyb3TQBgMKB/tNonXoxfnDLDEz+t6+Ox1P/+9hKjIxOAGD+RcvEOYYg+H5dMak/HwR25JwwGWd7amTu5FnPkZPjv9/1zgn0P3kfJqffHk+x60p9iM8PP+ToosXsCkMx4WGcvemHQYnJ9mC2CqKTzB5jx+ETJR+zp24VGuTc+M+jqQbc8X+2Onpem1IVKd7il34YlIgmTDiqKy6I1LrsMUqprsoUoGrNjsilBlVT2rF1lbPnSdLcgXV4Rw2Ak535EyaoBkRCKcSVpIyQUoqlNwbRSeYeI2v+jEsxf8Yl4xNpTV/Z7DnFzqZUxcRhCiIAjsjzJ0xEgKEzgBoJvBAXR3XFR99BTLAWHqNk6MPje2nK3EGkd3W2jHSZggiAI/L81CmlMp044OxD4aSU7thjdlStxJbUOkzGB4kd1RW7mAWRWhd2el42u2LvuqXYuvo6bSGxSO/qbBnpMgURAEfk+dwuVXZSbGjsMReqNBZWHURX3eN4dOp3rcsI8aRgVHdt09fw8lHkjb6DSK2La3peoUjv6mwZ6TIFEQBH5PncrlBz8mHPeUwNRtBZ16d1VBepglHdhtR6sxe2mKLEnV6ksXpbRrq69+81BEfkudzmpToZwVtabjcQBRe61MAebH1Q38g56swhz89X4k7PT6zedVs40rUK88j9cDIhZMukURgMy9GPOh/c8/P5LSsbZFvIKMwjD4OTEbzhq89CHaUaNqqLOnPI8/OFcBdX2JadvSfQtvZFrmmICcbIEy7UTAjD4pdR54N7fr4QNkzObQsAjCpoW9OQpHLJUWFH7lO5D6UNH1ht+e1OUjcDFnUVQc/PF8IFMLctVfLx73WsaeDiuOAxRu5TudijDXFJbW00LH6eJLo/l21rX8RHgyPjP19UV42965ZG9vw2KxUj54jcp3IjWhtWc2qrdW1LnnIUIr470V3fnCUPgsfJTp/KLTG3ob63tsU0MU3L9DR57HcXe5d0L6AypeRBnDC04lO5L67TL3UiKyPGNC3TU9gihHRDiiemH4ak3OjG6cgn9F3STWR4WqZXnsJpPu9OEjkQoDyMkRvAhlg6OeMp/usz3dBPFogNmVWuaciI0o0duQEKc3zPDo3E50sVBYO+uJ4mEn2mG/oZCLi+CPg511G9T7ZUbgwQO3IDZL/82fzeUQXm17phwBc3O7Jd8tjLAIBt990QWXE0P1kgri8Cfs51VO9TAjOi2JEbIBtLvyBVPf47hlhccPjFDTOMoHORi590QtcXAT+dZFQdrC2VGwPEjtwUp/uxNfVX2FO3Cptr16Gp6gMj0xWN5PCL67ezLXch0DnP4WfjC9cXAT+dZFQdbAglDkzH9ENTPLUU6sgOiEpjWFXhQGouGr/+MrMPnHCYyuh3RaHtq3gD4SdtNKYpp1FKVvqhjR+Y93ZDxm47a2UUrdILWNyJl0qJCyVVzmEqo98FWuVG3UEucjE6ndBP2mhMU05NEEhoRUSWisjbInJQRP4yiGP6YsDkl2sxi+uVCmPYGksGyseTg9zXk0WlyC3fI3IRqQbwBIDPAOgDsENE/kUptc/vsT0zcda60l2CptrdYY3+So1eTYglexXV0nKuKyC3ggitLAJwUCn1XwAgIv8I4GYA+jpyE+t4VKqnoem2M6xVpaXCGDbUnyklqholNp8j0iOI0Mo0AEdyfu4b+10eEblLRLpFpHtgYCCApy3DxFlrE+8SEN7or1QYQ3flPRvwHJFbkU12KqU2AtgIZLJWQn0yEydVTLxLQHijv1KjV92V94rxGl4KKyxl4jkiswUxIj8KYEbOz9PHfke5TLxLgBmjv7ILdZwu6/ax/NvN5GJuW294bBt29nJSkvTznUcuIjUAfg3g08h04DsA/F+l1N5S/03s8shtTHc0SG4OdpUAF6ZqoKDQOvVibKpeh9S73ZV3EvKx45Cb/PLcthbiTjcUttB2CFJKpQH8CYCfAdgPYEu5TjyWbEx3NEhunH5UAacH0+Oj3PS7bzibW/AxB3Ft0wi2pNZhT90qbEmtw7VNIyUfm9vWXJyUJJ0CySNXSv2bUupKpdQVSqm/DuKYVjF0ItMWhdUfs0ZGgX2jM53l1/vIw9+QWo+FVQfRIOewsOogNqTWO2prtQAXpjKbGSsFDI8ofRUrDaoASdFjrZUgxGwxT1Sy8eb33u3F5tpMnZnn6v4KTfIBgMwot6vpYWdzCz7mIFIDe1CDzCi8BiNIDewp+di8OYWZE3HlZRMgkrmT2H30ZLhlYMvhXWGixXOJftQ0LeaxXXaS8R+qv4tr5ABqZRQdchDPNHwfy4YeRuvUi/E3yzuACb9X9jiZ7JGD2HdsTSZ7ZFkHGhtcZI+4yCgqzChpW/uitzKwQe/P6eWukHM7scEReRB8bgyQVNl4c1tVL2ol0xuKSqNVel0tdfe9pD1nND/0iU780dmvOS51G2kZ2HK83BVyFB+pMMsosyMnbbKd4N7RmRhWYx9FD6Ep34uaci7Ey0fW4uWjcHxRiLQMbDleQkuc24lUmDV0khVa0X0rqfv5DZOtXfLnx/4MXan1uEp6UeUhNOV1UVOxBT2uLgqn+9G4dQW2/s9uYGb2/axwBxFWGM7LIjiHIaXc83Rt0wg2pNZn5hH4GXYlzBo6yapH7iPXOBbPH1NeV1gWqyEOwHld8Zz3M41q9KjZeHTqd8s+v1Elah0OLHLP05bUOiysOpiZHOZn2JUgatYnqx55KbpvJXU/f0x5XdJebIS07b4bnFc4zHk/azCCq3C4YuGxsIqUeeJwFJ97nlqldzzDh59hd8Ksnpmsjlx3vZMiz2/UCM2NGISJ8kIyAggESx57Ga1TL8a2+26o/D7kvJ/Dqgp7R2dWvGW2sURt7nnap2ZioeSMyJlq61iYNXSSNdmpu95Jkee3dhMBhxkPYc7U+5U7UXlBqgZnhtLu3oex9/OMXIDX1WzcM7ymYozez473oamQ2557njZethaj08yrGZR0yYqRG8jvPpJ++LobeGRaphPPSjVk0i8L2LKXpZ/3wc15NPIOTMfcTQzu6HRgjNxQOjcR8BWvdRimsiWU4Od9cHPLbGSJWh1zN2EtjEqoZIVWDKSzjKyvTtZhmMrIUEIRJpTz1UZHiQlO/AeKI3LNdI7QfN0NOMx4iGqfS7+MHClHRUeJCd2JBzHDGLkLRsY3fYjb6yGLMEbuSakYOTtyF2yZuCOieAptYwktNNVetmXijoiSxc6OXFPVNlsm7ogoWezsyDXNeCc6s4GIjGVn1oqmGe9EZzaEqNikKwC7JmIDmrxLwgR0El5j1Oyc7OSMt7G8fEn9ViE0omMIaHVkoBPqhn5PmDTgXbwmO7kjj7G81I4pNonsZmLZiHo1AYX7Ap1QN3QHICYNBM/OjpyM5eVLWmwS2c3EshEdQ0CrI8u9btcFyAxdPcmkgeCxI6dAefmSFptEdjOxbETHEFBlzXKv2/WdR1RL712mAzNpIHh2xsjJTKf7MfSTLyH97hvYNzoTXU0P429W3Bh6vNqIGHkEXFdojCpGzp2vIsPqhxS+LSuQercbKZVGZ9VB/PCCvwcm/F7oTxtaNpFhk4Wua+N42cfTi4BDOEm5MAeJoRUKTgQx2Ug3qjBsstDYkETAIRwjJq8twxE5BSeC/P5I97w0bLLQ2HUMAVdPNGLy2jLsyCk4EZRDjfRLzlKrzgQcwtG52Yqt2JEnWOCxyAhispF+yXXU6SZratibhFkrCWbSCjunF5WkTYQl7fVSecxaofOYEoscODWIGx7bhtODmXj0zt7SsW9j48QhiXROgKzFrJUEM2IhDTKdVbYTB4BRFcMJLo819E252JLZ2JF7EcbGFuWOGdJGGqaksxXrnAQSTYphVDymMppysSWz+YqRi8gyAN8GMBfAIqWUo8C39THyMFaylTtm3t+qgdSFgFJGLFIJwrKu7djVewIjYx/FKgEEwIjSH7sPzCPTMp14VqoBePBoxRg4Y+SUK6zqh3sA3AbgFZ/HsUsY+cXljpn3txFg8FTki1TCXIizYXkHFs7M3Bn8ZstE1NdWjXfqsQknlFg0U2nxS3ZOYO+6pdi6+jp24lSUr45cKbVfKfV2UI2xRhjFiModM/dvuSJcpBLmarvCzqrtE5fEL5xQoqgWY+AUBMbIvQio0p3jY+b+rW4CgOrM7yNcpBJlh2NK7D5QJWroMwZOQagYIxeRlwBMKfKnh5RSL4w9ZhuA+8rFyEXkLgB3AUBzc3NHb2+v1zYnm6ZCTiblnMcJY+DkRqkYeSALgpx05Lmsn+xMoLeOfYjf73oVpwfTaKirwXOrr8VVHD0SRSpeW71ROMqkOX7rhT04O5yZcD07nMa3Xtijq5XkUaSVIylSvjpyEblVRPoAXAvgX0XkZ8E0i7Qok+ts4qQcOyZ3WB42vvxmrfyzUmq6UqpOKXWZUupzQTWMNCiTAmnipJwpHZMtFxQTL8YUDIZW6GNlUiBNzCQxpWOK9ILiY5WviRfjwIW0Ctp07MjpY2VSIE1cmGJKxxTpBcXHrkUmXowDZ9iuTlFh9UP6WFR7PAbElLrVkdZI97GqOBGVIw3b1Skq7MjJWqZ0TIFcUJyuD+CuReUl9PxwYwkiEzgtxKZpQZjxsufl2BtAVTWgRoEpV8fu/HBjCSKTOQ0JWBb+ikw2Nh5kRVKLcLLTBgmdidcp8pTCMAqxJUlCY+NZ7MhtkNCZeJ0iz1EPoxBbkiT8QsjQihdRxykTPtrQIfIcdYZM/Ln92fO/kwnCjtyL3HhcdoQc5pcwoTPxOkWaUkj+JfxCyNCKF1GPkHnbHblELJ6hPLaUWiiG6YdehLFnZxiYqkbkmA0191nGNki2jJCTOknKLB/ywJTaPV4wRu6FLfG4pE6SRj2HQbFg87wIR+RxltSUrKRewErhHYojNs+LcEQeZ0lNyWKWTz7eoThiSu0eL9iR28Tt5KUtISAg2InZpF7ASuEdSuwxtGKTOE9eBvnashewB49m/k16pk5SQ2wJwo7cRKVimnEeWcX5telmS5YVecbQiolKxTTjHPt18NoGTg2eV/fbhJ2KnNLWfptCbOQJR+QmKjU6jfPIysFrM2WzZa9sbz+ZiyNyE5UancZ5ZOXgtdm8YAOwv/1kLo7I3YoiJ9fkkbfGnGRTNlv2yvb2k7lYa8UtW+qshEXj62eMnJLO/q3eTCkAFefsCifnWOPrt3nBBmB/+8lc9oRWTMmhjnNOrpNzHOfXT2QpezpyU0bCJsev/XJyjuP8+oksZU9oxZQc6jhnjjg5x3F+/USWsmdEzpFg+HiOiaxkz4icI8Hw8RwnT9RJBKYkLcSMPSNyIgpe1EkEpiQtxIw9I3Jyx+/IhyOneCp8X997M9okAlOSFoqwOc+fI/K48jvy4cgpngrfV6mKNp3U4PRVm2vhsCOPK78jH4NHTuRD4fs6OhLtBLfBE+o218LxFVoRkb8F8HkAQwAOAfhDpdTJANpFfvlN1zQl3dNUtoaeCt/XqddEO8Ft8IR6kjdf/jmAeUqpqwH8GsAD/ptEgfA78jF45ORKWEW+bA09xeV9DYHNmy8HVjRLRG4F8PtKqeWVHmt10Syyi4siX64mux6ZlunEs1INma3liEJUqmhWkDHyVQD+vUwD7hKRbhHpHhgYCPBpicpwEet3Ndll8KQdJU/FjlxEXhKRPUX+d3POYx4CkAawqdRxlFIblVKdSqnOxsbGYFpPVImLDtfVZBdDFHkGTg1iWdd2tK19Ecu6tmPg1KDuJiVKxclOpdSN5f4uIisB3ATg00pHcXOicm5/9vxJyRJcTXYZPGmnQ/ZuZmQU43czLNkbHb9ZK0sB/AWATymlzgTTJKIAuehwNyzvOC9GTs7YnLoXB35Xdn4fQB2An4sIALymlFrtu1VEGnDjB+9sTt2LA18duVLqN4JqCBHZi3czerHWChH5xrsZvbhEn4jIcuzIiYgsx46ciMhy7MiJiCzHyU4iw9m84QFFgyNyIsPZvOEBRYMdOZHhuGqSKmFHTmS41qkXo3rsm8pVk1QMO3Iiw9m84QFFg5OdRIbjqkmqhB05hY5ZF0ThYmiFQsesC6JwsSOn0DHrgihc7MgpdMy6IAoXO3IKHbMuiMLFyU4KHbMuiMLFETkRkeXYkRMRWY4dORGR5diRExFZjh05EZHl2JETEVlOlFLRP6nIAIBeBw+dDOD9kJtjE56P8/Gc5OP5yBe38zFTKdVY+EstHblTItKtlOrU3Q5T8Hycj+ckH89HvqScD4ZWiIgsx46ciMhypnfkG3U3wDA8H+fjOcnH85EvEefD6Bg5ERFVZvqInIiIKmBHTkRkOeM7chH5WxF5S0TeFJF/FpFLdbdJJxFZJiJ7RWRURGKfVlWKiCwVkbdF5KCI/KXu9ugmIk+JSL+I7NHdFhOIyAwReVlE9o19X/5Ud5vCZHxHDuDnAOYppa4G8GsAD2huj257ANwG4BXdDdFFRKoBPAHgdwC0AvgDEWnV2yrtngGwVHcjDJIGcK9SqhXAYgD3xPkzYnxHrpT6f0qp9NiPrwGYrrM9uiml9iul3tbdDs0WATiolPovpdQQgH8EcLPmNmmllHoFwP/obocplFLHlFK7xv7/KQD7AUzT26rwGN+RF1gF4N91N4K0mwbgSM7PfYjxl5T8EZEWAAsA/KfmpoTGiK3eROQlAFOK/OkhpdQLY495CJnbpU1Rtk0HJ+eDiCoTkQYAzwNYo5T6UHd7wmJER66UurHc30VkJYCbAHxaJSDxvdL5IBwFMCPn5+ljvyMaJyK1yHTim5RS/6S7PWEyPrQiIksB/AWALyilzuhuDxlhB4DZIjJLRFIA7gTwL5rbRAYREQHwJID9Sqm/092esBnfkQP4PoAJAH4uIj0i0qW7QTqJyK0i0gfgWgD/KiI/092mqI1Nfv8JgJ8hM4m1RSm1V2+r9BKRnwB4FcAcEekTka/obpNmvwXgSwB+e6zf6BGR39XdqLBwiT4RkeVsGJETEVEZ7MiJiCzHjpyIyHLsyImILMeOnIjIcuzIiYgsx46ciMhy/wvecg8BROTMLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "normal_pcs = []\n",
    "uniform_pcs = []\n",
    "for i in range(100):\n",
    "    normal_pcs.append(rd.normal(size = [100,2]))\n",
    "    uniform_pcs.append(rd.uniform(size = [100,2], low=-2, high=2))\n",
    "\n",
    "#we show the first point cloud of each type as an example\n",
    "fig, ax = plt.subplots()\n",
    "normal_pc = normal_pcs[0]\n",
    "uniform_pc = uniform_pcs[0]\n",
    "ax.scatter(normal_pc[:,0], normal_pc[:,1], s=15, label = 'normal')\n",
    "ax.scatter(uniform_pc[:,0], uniform_pc[:,1], s=15, label = 'uniform')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each of these point clouds, we are calculating the 0th-dimensional persistent homology of its rips complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgms_normal = []\n",
    "for pc in normal_pcs:\n",
    "    rips = gudhi.RipsComplex(points=pc).create_simplex_tree(max_dimension=1)\n",
    "    rips.compute_persistence()\n",
    "    dgms_normal.append(rips.persistence_intervals_in_dimension(0))\n",
    "    \n",
    "dgms_uniform = []\n",
    "for pc in uniform_pcs:\n",
    "    rips = gudhi.RipsComplex(points=pc).create_simplex_tree(max_dimension=1)\n",
    "    rips.compute_persistence()\n",
    "    dgms_uniform.append(rips.persistence_intervals_in_dimension(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see two examples, the barcode of a normal point cloud and of a uniform point cloud"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:title={'center':'Persistence barcode'}>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEICAYAAACtXxSQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQZElEQVR4nO3dS3IbyXbG8e84POtJXcqayhGlSY/R1AoM7QDSXYGpCC+AUq/gmtwBsQGHTA08J3bAx9wDIcJj3aAw8PzcQSXEUgnEo1B41Mn/L4LRQFUByCTELw+ystDm7gIA9Ns/HboBAIDtEeYAEABhDgABEOYAEABhDgABEOYAEABhDkmSmZVmdmFmbmb3Znaefq7M7NrMBjt63Sszu9rFc++bmQ3T7+760G15jpkV6f1kTXIwxjpz1JnZV0lX7n5Z21ZK+irptbtPO369oSS5+2SNY8/r7TpGZnYm6a27vzt0W5YxM3d3O3Q70B0qcywyq99JAf4g6WPXL+Tuk3WCPHnT9evvwOOhG4A8EeboBTO7OHQbgGP2z4duAI5fmmYZSPr3dL+Q9KekW1XV8o27T9KUyYWkO0k3kv4q6bO7fzGzc1XVfSHpbTquSP+Vu79Nz73ouDL9DNL+mbuPV7RlkB47kzSfk38r6ev8sQseP5NUuPuXZc+9xq+sMLNRun2S+ld/zUHaXjSfd8XvcFVbzyTNp8HKxlTZhaqpskc1Pnlt2VccC3fnh58fP6r+4K8ljVSFw3m6P2gcU9Tu38/vp8d8VRVUg/RzJmlUO340fz5JQ1XBoRXHjSRdP9Pe59oySvvLdL+o/sn/9Pj7BX0brHruJb+/kaoBo77tStJZ43nq/fzeOP6X3+Eabb1vtHWk6tyHVA0K9ceVC34PG/eVn+P6oTLHIreeKr6mVHHO3H1W23wn6b2ksdKccdr/kB5zIukq/Xfy3HOrqirXOW7dtsxSW6bzNpmZzKxIt4eqgv6h9vh37v6wxnMv89C4f51+5o97540TyfM2pbuLfoer2qp6W72q5K/N7L+aj3P3qdnTuc8t+4ojQZhjU6X0tAoludbTx3upEWZeTXt8kPRB0oWZTSX9WyM81j5uw7YsW31TNvfXQm+d517XVFWVPfc4ny6qPd+Jfp7+aA4Iq9q66MTrTNXU0mzBvuZzd9VXHAhhjk1NVX38Xns+1cyG6fj5vPCVqqmEy5bHlW3bsqAv5ZJ92zx3XaGfg/FeqaqWpHmV3KjON23ryTOveyvpP1a0r8u+4kBYzYKNpKmPWS1Q5xccDZc8bNDY/9xFQsuOq1e3pbtP12zLopCb92UiaVo/3swGZjZo2c+5Zuh+UDrROz/5WQvyonbcs8+9RluLRltHkr64+/+kxw3qj2s89zZ9xZHgoiFI+lHtjlSFzoOqFRQLL9BprHyQqvnW+WqWj6rC7ErSOM1Nn6Xj5lMBpbtfplD5U1WI1dew/3Rc7XWvVFW18sWrWeptma9mOU3P/d/puHNV88AXae64SMfdp9ed+dPKkoXPveTXOJ+qeNRToJeS5gNPvR9f9TSVUqiaDrlSNfj88justWedtp6oqrQvF+ybpde8Sb+Hj+k92rivOC6EOQAEwDQLAARAmANAACvD3MxGZnazYNuwNhe6cBsAYD9Whnnzwo3aBQo/Lj9etK37pgIAntNmmuWNntbMTlVdbrxoGwBgT9pcNFQ07r94ZttP0vTLmST99ttvf/z+++8tXhoA8nV/f/93d3+5aF+bMJ/p1wsxFm37SVoXPJak09NTv7u7a/HS7Xw///Tsvr9c/ufe2gEA2zCz/3tuX5swv1XtSjxVFx8UC7YBAPZkZZink5mnZjZy9y/p29jO0/aidtLzl2198FzVTsUOoE9WhnkK5r80ts0vsZ4s2wYA2A8uGnrG9/NPS+faAeCYEOYAEABhvgLVOYA+IMwBIIAs/k9DrEwBEF0WYb7tVAmDAYBjl0WYE8YAossizLs4icmAAOCYZRHmBDGA6LII810sL2SAAHBMWJoIAAEQ5gAQAGEOAAEQ5i3xRVwAjglhDgABEOYAEABhDgABZLHOnDXhAKLLIsx3eaKSgQLAMcgizAlcANFlEeZU5gCiyyLMCVwA0WUR5ru+uIfBAsChsTQRAAIgzAEgAMIcAAIgzDvAF24BODTCHAACIMw7wlfiAjgkwhwAAshinTnrwAFEl0WYH2r6g0EEwL5kEeaEKoDosgjzQ56YZCABsA+cAAWAAFqFuZmNzGxoZmfLtoELigDsx8ZhbmZDSVN3n0iamtnAzEaSlLbNjwEA7EmbOfM7Sfdm9k5S6e4TM/urpM9p/1TSQNKkozb2Xr06Zw4dwC5sXJm7+0zSlaRrSX+kzUXjsBfNx5nZmZndmdndt2/fNn1ZAMASbadZJu7+WtIsTbHMJJ0se5y7j9391N1PX7582aqxETCHDmAX2pwAHbj7Q7r9N1Uhfqun6ryUdLN90wAA62oT5uM0ZTKU9D5V3F8klWlbMT8RisWozgF0beMToGnOfLxg+2W6SZADwJ5lcQUoK0gARJdFmB/jtAYDDIAuZRHmBCeA6LII82OszCUGGQDd4Yu2ACAAwvyAjvUTA4D+IcwBIADC/MCozgF0gTAHgAAI8yPw/fwTFTqArRDmABBAFuvMWc8NILoswrzPUxgMRADWkUWYE4gAossizKnMAUSXRZgTiACiyyLM+1yZSwxGAFZjaSIABECY90DfP1kA2D3CHAACIMx7guocwDKEOQAEQJgDQACEOQAEkMU6c9ZpA4guizCPdvKQwQlAUxZhTvgBiC6LMI9Wmc8xSAGYyyLMCT0A0WUR5lTmAKJjaSIABECY91jUTxwANkeYA0AAhDkABECYA0AAhDkABNBqaaKZDSSVkuTuX9K2kaSZpNLdx101EACwWtt15n+6+zszOzOzUtJAktx9krYN3X3SXTO3w3psANFtHOZmdibp1sx+VOBm9kHS53TIVFW4TxY87kySXr16tU2bNxZ9CR+DFYA2lfnr9N9HM7uS9FFS0TjmRfNBKfjHknR6euotXrc1wg5AdG2nWb66+8zM7lVV2zNJJ521qmPRK3OJAQvIXZswv9VTcBeqgnyqp+q8lHSzZbsAABvYeGliWr1SmNkw3R+nbWXaVhzTyU8AyEGraRZ3v0w3J8u2AQD2g4uGgvh+/imLcwMAFiPMASAAwjwYKnQgT4Q5AARAmAdFdQ7khTAHgACy+B86c3UkgOiyCPNcpxwYxIB8ZBHmhBqA6LII81wrc4mBDMgFJ0ABIADCPLicP5UAOSHMASAAwjwDVOdAfIQ5AARAmANAAIQ5AASQxTpz1loDiC6LMOcE4M8Y3IB4sghzwgtAdFmEOZX5YgxyQBxZhDmhBSC6LMKcyvx5DHRADCxNBIAACHMACIAwB4AACPPMcT4BiIEwB4AACHPo+/knKnSg5whzAAggi3XmrKUGEF0WYc4UwvoY+IB+yiLMCSgA0WUR5lTm62PgA/opizAnoABEt1WYm9mFu39Mt0eSZpJKdx930LbOUJlvhsEP6J/WSxPNbCipTLdHkuTuk9o+AMCetApzMyslTWub3tTuTyUNtmwXDoiLiID+aVuZl+5eD/Oisf9F8wFmdmZmd2Z29+3bt5YvCwBYZOMwN7PhfDqlZibpZNnj3H3s7qfufvry5ctNXxYAsESbE6CPaU68kFSa2UDSrZ6q81LSTSetAwCsZePK3N0fUmV+ohTg7v5FVbAPJRULKnf0EPPmQH+0XpqYlh+Oa/cv002CHAD2jIuGACCALMKc6YLdYJAEjkcWYU7oAIguizCnMt8dBkrgOPB/GgKAAAhzbIVL/4HjQJgDQABZzJlj9+rVOfPowP5RmQNAAIQ5Osc8OrB/hDkABECYY2eo0IH9IcwBIIAsVrOwugJAdFmEOR/1D48BFditLMKcIAEQXRZhTmV+eAyowG5xAhQAAiDMsRcsUwR2izAHgAAIc+wVFTqwG4Q5AARAmANAAIQ5AASQxTpz1jgDiC6LMOeE23FikAW6k0WYExoAossizKnM+4XBF9hcFmFOOACILoswpzLvJwZhYH0sTQSAAAhzAAiAMAeAAAhzHC2+lAtYH2EOAAFksZoF/fZcdc5qF+DJxmFuZoWkMv28cfePaftI0kxS6e7jDtsIAFihTWX+XpLcfWxmb8zsTNJj2jYxszMzG7r7pMuGboMKDkB0G4d5o+ouJd1Ieivpc9o2lTSQ9FOYp9A/k6RXr161aWtrnESLjwEbuWs9Z25mpaTHVI2/a+x+0Tw+DQJjSTo9PfW2r9sGf+gAotvmBOjI3T+k2zNJJ9s3ZzeozCExqCO2VksTzWzk7pfp9kDSraQi7Z5PvQAA9mTjMDezoaQLM7s3s3tJJ+7+RVKZ9hXHdPITmOMTGiJrcwJ0Iun1gu2X6SZBDgB7xhWgABAAYQ4AARDmyArz5oiKMAeAAPiiLWSnXp2z9hxRUJkDQABZVOZUXwCiyyLMOemFZRjsEUEWYc4fK4DosghzKnNsgsEffcQJUAAIIIvKHNgE/89R9BGVOQAEQJgDa+LcC44ZYQ4AARDmwAa+n3+iQsdRIswBIIAsVrOwCgFAdFmEOR+LsU8UDziELMKcPy4A0WUR5lTm6AsKD7SVRZjzBwIguizCnMockVGsQGJpIgCEkEVlDkS26pMnlXseqMwBIAAqcyC4Q54z4lPB/lCZA0AAhDkABECYA0AAWcyZM28HILoswpyLhoDDo6jarSzCnH9EAKLLIsypzIH85FbEZRHmub2pAPLTWZib2UjSTFLp7uOunrcLVOYAduGYCsVOliamIJe7T9L9YRfPCwBYT1eV+RtJn9PtqaSBpElHzw0AR+mYPvV3FeZF4/6L5gFmdibpLN39fzP7345eex3/Iunve3y9Q8upvzn1Vcqrvzn1VVqvv//63I6uwnwm6WTZAWke/SBz6WZ25+6nh3jtQ8ipvzn1Vcqrvzn1Vdq+v11dzn+rp+q8lHTT0fMCANbQSZi7+xdJZTrxWcxPhAIA9qOzpYnufpluHmOQH9VSyT3Iqb859VXKq7859VXasr/m7l01BABwIHwFLgAEQJgDQADhvptl1dcKHPPXDmxqWV/MrFC1sqiU9MbdP+69gR1b970zs4u+93eNf8cDVe/tfAFCr+X0dyv96M8Hd3+7ZP9MG/Q3VGW+6msFIn3twBp9eS/pdP6Hni7a6q1137u0vdxj0zq3Zl//TO/tiZmF7m+6P037p2kg67VlA3DbnAoV5qq+VmCabs+/VmCT/X2ytC/uPq6N6GXt2L5a+d6lUOt7P6UVfU0D862Zlel97nufV723d5Ku559G3P1hn407gFY5FS3Mi8b95tcKrNrfJ0Xj/sK+pIB7DLD2v2jcX9TfMkCwSav7+jptezSzqzSl1mdF4/5P/XX3maQrSdeS/thPkw6qaNxfK6eihflMy79WYNX+Pplpvb6M3P3DjtuyDzMt6a+ZDQMMWHMzrX5vv6aQu9fTdx711Uwr3ltJE3d/LWk2n4YIbKYWORUtzFd9rUCkrx1Y2RczG80v5gowz7iqv49mNkx/6GXP+7vOv+O5QtUff5+t6u+gNrXyN8UpyJ7TKqdChflzXytgZjfL9vfRqr6m7Rdmdm9m9+r5H8Aa7+1D2naiXz+m9sqa/46L+Ymxvq/uWNVfSWMzO0v73/e9v9KPv8/T+qeMbXOKK0ABIIBQlTkA5IowB4AACHMACIAwB4AACHMACIAwB4AACHMACOAfTCignQ5iR7UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gd.plot_persistence_barcode(dgms_normal[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:title={'center':'Persistence barcode'}>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEICAYAAACtXxSQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARHklEQVR4nO3dTXLb2HrG8edNZdYTXDk9daroSY/Z9ApC74DurODSVVmA7F5BIu1A2kCqrzzIXNqBPuYZmFUZ+5bMQebvHeCwBcP8AEAQH+f8f1UskyBIvqTIBy8PDmhzdwEAxu2f+i4AAHA8whwAIkCYA0AECHMAiABhDgARIMwBIAKEOSRJZjYxswszczN7NLPzcLoysxszm57oca/M7OoU9901M5uH1+6m71p2MbMs/D2ZkxwZY545iszsi6Qrd78sLJtI+iLpjbuvWn68uSS5+12Fdc+LdQ2RmS0lvXP3933Xso+Zubtb33WgPXTm2GZdvBAC/EnSx7YfyN3vqgR58Lbtxz+B574LQJoIc4yCmV30XQMwZP/cdwEYvjDMMpX013A5k/S7pHvl3fKtu9+FIZMLSQ+SbiX9u6Q/3P2zmZ0r7+4zSe/Celn4V+7+Ltz3tvUm4TQN16/d/fpALdNw27WkzZj8O0lfNrfdcvu1pMzdP++77wovWWZmi3D+LDy/4mNOw/KsfL8HXsNDtS4lbYbBJqWhsgvlQ2XPKn3zOvK5YijcnROnP0/KP/A3khbKw+E8XJ6W1skKlx83l8NtvigPqmk4LSUtCusvNvcnaa48OHRgvYWkmx317qplEa6fhMtZ/pb/7vaPW57b9NB973n9Fso3GMVlV5KWpfspPs9vpfV/eA0r1PpYqnWhfN+HlG8UirebbHkdaj9XTsM60Zljm3sPHV9Z6DjX7r4uLH6Q9Juka4Ux43D9U7jNmaSr8O/drvtW3lVWWa9qLetQy2pTk5nJzLJwfq486J8Kt3/v7k8V7nufp9Llm3Da3O69l3Ykb2oKF7e9hodqVbFWzzv5GzP77/Lt3H1l9rLv88jnioEgzFHXRHqZhRLc6OXrvVQKM8+HPT5I+iDpwsxWkv6tFB6V16tZy77ZN5Py9YXQq3LfVa2Ud9kbz5vhosL9nen74Y/yBuFQrdt2vK6VDy2tt1xXvu+2nit6QpijrpXyr9+Vx1PNbB7W34wLXykfSrhsuN6kaS1bnstkz3XH3HdRpu+D8VGhq5akTZdc6s7r1nq243HvJf3HgfrafK7oCbNZUEsY+lgXAnVzwNF8z82mpet3HSS0b71idztx91XFWraF3Oa53ElaFdc3s6mZTRs+z41y6H5Q2NG72flZCPKssN7O+65Qa1aqdSHps7v/T7jdtHi70n0f81wxEBw0BEl/drsL5aHzpHwGxdYDdEozH6R8vHUzm+Wj8jC7knQdxqaXYb3NUMDE3S9DqPyuPMSKc9i/W6/wuFfKu1r59tksxVo2s1lm4b7/FtY7Vz4OfBHGjrOw3mN43LW/zCzZet97XsbNUMWzXgJ9Immz4Sk+jy96GUrJlA+HXCnf+PzwGhbqqVLrmfJO+3LLdevwmLfhdfgY/ka1nyuGhTAHgAgwzAIAESDMASACB8PczBZmdrtl2bwwFrp1GQCgGwfDvHzgRuEAhT8PP962rP1SAQC7NBlmeauXObMr5Ycbb1sGAOhIk4OGstLlVzuWfScMvywl6aeffvr1l19+afDQAJCux8fHv7v7z9uuaxLma/14IMa2Zd8J84KvJWk2m/nDw0ODh27m2/mng+v85fK/OqgEAJozs//bdV2TML9X4Ug85QcfZFuWAQA6cjDMw87MmZkt3P1z+DW287A8K+z0/GHZmFTp3svo5gEMxcEwD8H8l9KyzSHWd/uWAQC6wUFDR/h2/qlRRw8AbSPMASAC/J55C7Z154ynA+gSnTkARCCJzpwuGUDskgjzvndSsjEBcGpJhDlhCiB2SYR53515ERsWAKeQRJgToABil0SYD6kz34aNDYBjMTURACKQRGc+dMVvDnTpAJqgMweACBDmABABwhwAIkCYDww/qwugCcIcACKQxGwWZogAiF0SYR7zsAUbKgBSImFO4AGIXRJhHnNnXsaGC0hTEmFOwAGIXRJhnlJnLrHxAlLE1EQAiABhDgARIMwBIAKEOQBEgDAHgAgQ5hHix7qA9BDmABCBJOaZM+8aQOySCPOUhxzYkAFpSCLMCTQAsUsizFPuzCU2ZkAK2AEKABFo1Jmb2ULSWtLE3a93LcMw7PpmQscOxKN2Z25mc0krd7+TtDKzaQhyhWWbdQAAHWkyzPIg6cbMpsq78CdJbyWtwvUrSdOW6sMJcXAREI/aYe7ua0lXkm4k/RoWZ6XVXpVvZ2ZLM3sws4evX7/WfVgAwB5Nh1nu3P2NpHVhrPxs3+3c/drdZ+4++/nnnxsVi9OgQwfGr8kwyzQMrUjSfyoP8Xu9dOcTSbfHlwYAqKrJbJZrM1sqHxsvzmY5D117ttkRCgDoRu0wD2PmP0w9dPfLcJYgB4COJXEEKPOpAcQuiTBn5141bPSA8UoizAkpALFLIszpzKtjwweMEz+0BQARIMwBIAKEOQBEgDDHdzi0HxgnwhwAIkCYYyu6c2BcCHMAiEAS88yZOw0gdkmEOUMGzbARBMYjiTAnlADELokwpzM/HhtEYNiSCHOCCEDskghzOvPjsUEEho2piQAQAcIcACJAmANABJIYM8fxyvsdGEMHhoXOHAAiQGeORoqdOl060D86cwCIQBKdOZ0jgNglEeYcNNQNNppAf5IIc0IGQOySCHM68+6w4QT6kUSYEzAAYpdEmNOZd4uNJ9A9piYCQAQIc7SOb0JA9whzAIgAYY6T+Hb+iQ4d6BBhDgARIMxxUnTnQDcaTU00s6mkiSS5++ewbCFpLWni7tdtFQgAOKzpPPPf3f29mS3NbCJpKknufheWzd39rr0yj8O8ZwCxqx3mZraUdG9mf3bgZvZB0h9hlZXycL/bcrulJL1+/fqYmmvjq/5wsGEFTqNJZ/4m/PtsZleSPkrKSuu8Kt8oBP+1JM1mM2/wuI0RIABi13SY5Yu7r83sUXm3vZZ01lpVLaMzHxY2rkD7moT5vV6CO1Me5Cu9dOcTSbdH1gUAqKH21MQweyUzs3m4fB2WTcKybEg7PwEgBY2GWdz9Mpy927cMANANDhpC59iHAbSPMAeACBDm6AU/xAW0izAHgAgk8d/GYbi2defMQwfqozMHgAgk0ZnT6QGIXRJhzo628WEDDNSTRJgTDABil0SY05mPDxtgoB52gAJABAhzDBIHFQH1EOYAEAHCHINGdw5UQ5gDQASSmM2CcdvVnTPjBXhBZw4AEUiiM6eDAxC7JMKcnWhxY2MNJBLmfNgBxC6JMKczjx8bbKQuiTDngw4gdkmEOZ15Wth4I0VMTQSACBDmiA7fxJAiwhwAIkCYI0p050gNYQ4AESDMASAChDkARCCJeebMOwYQuyTCnJ1h6WJDjlQkEeZ8oAHELokwpzOHxEYdcUsizPkQA4jdUWFuZhfu/jGcX0haS5q4+3ULtbWGzhxFbNwRo8ZTE81sLmkSzi8kyd3vCtcBADrSKMzNbCJpVVj0tnB5JWl6ZF0AgBqaduYTdy+GeVa6/lX5Bma2NLMHM3v4+vVrw4cFAGxTe8zczOab4ZSCtaSzfbcL4+jXkjSbzbzu4wJtKe5DYfwcsWiyA/Q5jIlnkiZmNpV0r5fufCLptpXqAACV1B5mcfen0JmfKQS4u39WHuxzSdmWzh0AcEKNpyYWh03C5ctwliAHgI5x0BAARCCJMOegIezDxh4xSCLM+bACiF0SYU5njn3Y2CMG/E9DABABwhzJ45sbYkCYA0AECHMAiABhDgARIMwBMW6O8SPMASACScwzB6ood+fMP8eY0JkDQASS6MzpsADELokwZ+cWmqIRwFgkEeZ8IAHELokwpzNHUzQCGAt2gAJABAhzYI9v55/4ZodRIMwBIAKEOVAB3TmGjjAHgAgQ5kBFdOcYMsIcACKQxDxz5goDiF0SYc7XY7SJ5gBDlESY8+EDELskwpzOHKdAk4AhSSLM+dABiF0SYU5njlOgScCQMDURACKQRGcOnELxGx9dOvpGZw4AESDMgRawXwZ9I8wBIAKEOdASunP0qfYOUDPLJE3C6a27fwzLF5LWkibuft1ijQCAA5rMZvlNktz92szemtlS0nNYdmdmSzObu/tdm4Ueg5kGAGJXO8xLXfdE0q2kd5L+CMtWkqaSvgvzEPpLSXr9+nWTWhvj6y+6RgOBrjWeZ25mE0nPoRt/X7r6VXn9sBG4lqTZbOZNH7cJPlgAYnfMQUMLd/8Qzq8lnR1fzmnQmaNrNBDoWqPZLGa2cPfLcH4q6V5SFq7eDL0AADpSO8zNbC7pwswezexR0pm7f5Y0CddlQ9r5CfSBb4PoWpMdoHeS3mxZfhnOEuQA0DEOGgJO5Nv5Jzp0dIYwB4AIEObAidGdowuEOQBEgDAHOkB3jlMjzAEgAkn8t3EcjQcgdkmEOV9xMTY0IKgriTDngwEgdkmEOZ05xo6GBIewAxQAIkCYAyPAt0scQpgDQASSGDMHYrCrO2c8HRKdOQBEgTAHgAgQ5gAQgSTGzBlTBBC7JMKcaV2IHQ0Lkghz3ugAYpdEmNOZA7vR7MQhiTDnzQogdkmEOZ05UB3NzzgxNREAIpBEZw6guirfZOneh4fOHAAiQJgDqO3b+Sf2RQ0MYQ4AESDMATRGdz4chDkARCCJ2SzseQcQuyTCnK+CQP9oqk4riTDnTQQgdkmEOZ05gDrG2AAmEeZj/MMAQB2thbmZLSStJU3c/bqt+20DnTmAUxhSo9jK1MQQ5HL3u3B53sb9AgCqaaszfyvpj3B+JWkq6a6l+waAQRrSt/62wjwrXX5VXsHMlpKW4eL/m9n/tvTYVfyLpL93+HhNjaVOaTy1jqVOaTy1jqVOaTy1Vq3zX3dd0VaYryWd7VshjKP3MpZuZg/uPuvjsesYS53SeGodS53SeGodS53SeGpto862Due/10t3PpF029L9AgAqaCXM3f2zpEnY8ZltdoQCALrR2tREd78MZ4cY5IOaKrnHWOqUxlPrWOqUxlPrWOqUxlPr0XWau7dRCACgR/wELgBEgDAHgAhEF+ZmtjCzeZjXXvv6rlSscxCzgvbVamaZmU3DOhd91Feo5dBrOg+nQddZWK/XOkMNh17Tb2Z2a2bnXddWquNQnZv36KLr2rbUsu/zNDWzL2b2GE6V3wNRhfmhnxUYys8OVKkjzBDqXYVaf5M029Tb10aywt9+Kmkarp+a2aT7Kqu/B8PyXmos1FCl1vfu/q4wAaJzFev8PbxHz/r620uVaj1z9zfu/qukv0q6qnrfUYW58p8VWIXzm58VqHN9V4ZSRxV7a3X368IPq00K63btUJ1P7n5pZpmklbsPsk5JCmHTV31FVd6nWZ/hGOytMzQY92Y2Ce/XPl/bQ+/T4mzASZ1aYwvzrHS5/LMCh67vSla63FcdVWSly1trDR/o5x6PMchKl3e9pjNJX05byl5Z6fK2Omt9iE8oK13eVuuZpGczq9xBnkBWulyu801Y9mxmV2GD3pfyY+/6PC3rfjuPLczX2v+zAoeu78paw6ijirWq1bpw9w8nrmWftSrUGTY2b3ocO11rT51mNh/QQXdrVfiZDndfS1oP9TUNvoQ6H/XyG1F9WKva5+ld3TuOLcwP/azAUH52YCh1VHGwVjNbbMZMw9h0H/bWaWYXhfH8tfrbmB56PZ/DzrGF8qOq+xyCO/SaLnuub6PK534jU/7370uVz1NWXlZFVGG+62cFNrNChvKzA4fqDOfnkmZ9730/VGtYfrHZ+66eQrLCa3olaVW4vpcjAyu8R5/CsjP9+JW8UxVe07+Fy4vC+oOrM1yfbXY29vmf51T57CsMXdW9b44ABYAIRNWZA0CqCHMAiABhDgARIMwBIAKEOQBEgDAHgAgQ5gAQgX8Az37XSTlAh74AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gd.plot_persistence_barcode(dgms_uniform[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we can see that there are actually differences between the two barcodes. The length of bars from the normal point cloud seems to have a greater variance. We are using PE to proving this fact."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that in order to apply entropy, all bars in the barcode must be finite. A common way of proceed is to change infinity by a fixed value. Note that infinity bars does not contain information in the 0-th dimensional persistent homology of rips complexes, since every barcode has only one starting at 0. Then, in this case, we can remove the infinity bars without losing information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we create a lambda function which removes the infinity bars from a barcode.\n",
    "remove_infinity = lambda barcode : np.array([bars for bars in barcode if bars[1]!= np.inf])\n",
    "# apply this operator to all barcodes.\n",
    "dgms_normal = list(map(remove_infinity, dgms_normal))\n",
    "dgms_uniform = list(map(remove_infinity, dgms_uniform))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once all bars in our barcodes are finite, we can apply PE. We use the [representations.Entropy](http://gudhi.gforge.inria.fr/python/latest/representations.html#gudhi.representations.vector_methods.Entropy) method. When applied to a list of barcodes, it gives another list with their respective persistent entropys."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[29.29791597],\n",
       "       [29.51345627],\n",
       "       [29.25127001],\n",
       "       [29.67212005],\n",
       "       [27.79781263],\n",
       "       [30.88106544],\n",
       "       [25.05627269],\n",
       "       [30.77649594],\n",
       "       [22.95157941]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PE = gd.representations.Entropy()\n",
    "pe_normal = PE.fit_transform(dgms_normal)\n",
    "pe_normal[:9]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[31.59298874],\n",
       "       [29.9528389 ],\n",
       "       [31.21758435],\n",
       "       [30.24050121],\n",
       "       [32.23988175],\n",
       "       [30.86782984],\n",
       "       [30.66243524],\n",
       "       [31.41764578],\n",
       "       [31.01308312]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pe_uniform = PE.fit_transform(dgms_uniform)\n",
    "pe_uniform[:9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's show the differences using [boxplots](https://en.wikipedia.org/wiki/Box_plot). In order to do that, we need to transform these two ndarrays into two arrays of numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe_normal_array = np.array(pe_normal[:,0])\n",
    "pe_uniform_array = np.array(pe_uniform[:,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the boxplot function from the matplotlib library we obtain the following graphic,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD3CAYAAAA9vL6wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAM2klEQVR4nO3dsWtbZxfH8d953UKh0Cp6q7VvUCaBJlXZvASc/0BJ8CxwJq81RUPTwRR7K5ksyCxC8h9EoEVbZE8BT1FoV5fEFAoF05530KNEcWXJdnx1fezvB0Tvfa6ke0KVX46eR7oydxcAII7/5F0AAOBsCG4ACIbgBoBgCG4ACIbgBoBgCG4ACOazrE/wzTff+M2bN7M+DQBcKbu7u7+7e2nascyD++bNmxoMBlmfBgCuFDP79aRjTJUAQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDCZf44bwNViZud6HNf+vzgEN4AzOSmAzYxwXhCmSgAgGIIbAIIhuAEgGIIbAIIhuAEgGIIbAIIhuAEgGIIbAIKZ+wUcM1tJm3fdfcPMCpLK6Xbb3TcyrA8AcMzMjtvMapJq7t6VVDOzsqT7kuru/jzdZy37MgEAYzM7bnffk7SXuuyhuw8ltSfuUpb0IrvyAADHnXaOuy7p9eRA6r7fpm5cx46tmdnAzAYHBwcXUCYAYOxUwZ3C+ZaZNSaGG+7+8IT7t9297u71Umnqr8sDAM5p3hz31sQc9qGkYhpvuPt22q5lWiEA4CPzOu4dScP0yZKCu7fT9paZ7ZrZrlKYAwAWY97i5FDSMO1201hX0q2M6wIAnIAv4ABAMAQ3AARDcANAMAQ3AARDcANAMAQ3AARDcANAMAQ3AARDcANAMAQ3AAQz9xdwkA8zO/Nj3D2DSgBcNgT3JXVSCJsZAQ1cc0yVAJiqWCzKzE59k3Sm+5uZikUuLnoedNwApnr37l3m7+7OMyUIOm4ACIfgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBguKwrgKn8x6+kR19nfw6cGcENYCr76Y+FXI/bH2V6iiuJqRIACIbgBoBg5k6VmNlK2rzr7htprCHpUFLZ3dvZlQcAOG5mx21mNUk1d+9KqplZOYW20thksAMAFmBmcLv7nrtvm1lB0tDdh5JuSxqmuwwl1bItEQAw6bRz3HVJr9N24dix/x6/s5mtmdnAzAYHBwefUB4A4LhTBXeaFrk1MbddnHP/trvX3b1eKpU+vcorqlgsyszOdJN0pvsXizP/VwEIaObipJltSXqdFiAPNQrsl/rQdZclvciwvivt3bt3C/mcLICrZV7HvSNpmBYgC6mTfi6pPDHWzbxKAMB7MzvutBg5XojsToxvHx8DACwGX8ABgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhh8LBnCirC9SduPGjUyf/6oiuAFMddYrV5pZ5le7xAhTJQAQDB13jvzHr6RHX2d/DgBXCsGdI/vpj4X8kII/yvQUABaMqRIACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACGbmL+CYWUFSOd1uu/tGGm9IOpRUdvd2xjUCACbM67jvS6q7+3NJMrM1M1uRNHT3rqShmdWyLhIA8MHMjvtYN12W9ELSQNKumd3TqOPuZlgfAOCYU81xm1lZ0lt377r7oaQdSc8kfXfC/dfMbGBmg4ODgwsrFgBw+sXJhrs/lKQ0VdJ191uSDtN890fcve3udXevl0qlCywXADA3uM2s4e7babsmqebue+nwz5KKGdYHADhm3qdKViRtmdkPaWhDUtvM1iQNxadKgGvHzM51zN2zKOdamrc42ZV0a8ohwhq4pgjg/PEFHAAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIhuAGgGAIbgAIZua1SpC9WRfluQg3btzI9PkBLB7BnaPzXKzHzLjID3DNMVUCAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AAQDMENAMEQ3AA+SafTUbVa1dLSkqrVqjqdTt4lXXn8dBmAc+t0Omq1Wnry5ImWl5fV7/fVbDYlSaurqzlXd3XZrN8vNLOCpHK63Xb3jTReS2Ny9+ezTlCv130wGFxUvdcevzmJy6Rarerx48e6c+fO+7Fer6f19XW9evUqx8riM7Ndd69POzZvquS+pPo4nM1sLY3/kMaKZla+uFIBRLK/v6/l5eWPxpaXl7W/v59TRdfDzOB297a7t9NuWdIwhfdLMyun48PjjzOzNTMbmNng4OAgg7IBXAaVSkX9fv+jsX6/r0qlklNF18OpFidTV/3W3buSbkn6r6S3ZraTplM+kgK97u71Uql0oQUDuDxarZaazaZ6vZ6Ojo7U6/XUbDbVarXyLu1KO+3iZMPdH07sv3b3QzPblbQmafviSwNw2Y0XINfX17W/v69KpaLNzU0WJjM2N7jNrOHu22m7JumlpGI6XJB0mFVxAC6/1dVVgnrBZk6VmNmKpC0z203ddTEtShbSMU3MgQMAFmBmxz0xp318fDw10s2iKADAyfjmJAAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEQ3ADQDAENwAEMzO4zaxgZjUza5jZ1pTj/xoDAGRrXsd9X1Ld3Z9LkpmtjQ+Y2Yqkcoa1AQCm+GzWQXdvT+yWJb2QJDMrSxpmWBcA4ASnmuNOQf3W3btpqOzuJwa3ma2Z2cDMBgcHBxdRJwAgOe3iZMPdH0qjKZKJAJ/K3dvuXnf3eqlU+uQiAQAfzJwqkSQza7j7dtquSXqb5rcLkspmVnP3vWzLBACMzftUyYqkLTPbNbNdSUV330sdd1Gj8AYALNC8xcmupFsnHGtLak87BgDIDl/AAYBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBgCG4ACIbgBoBg5l4dEPkwszMfc/esygFwiRDclxQhDOAkTJUE0el0VK1WtbS0pGq1qk6nk3dJAHJCxx1Ap9NRq9XSkydPtLy8rH6/r2azKUlaXV3NuToAi2ZZvyWv1+s+GAwyPcdVV61W9fjxY925c+f9WK/X0/r6ul69epVjZQCyYma77l6feozgvvyWlpb0119/6fPPP38/dnR0pC+++EJ///13jpUByMqs4GaOO4BKpaJ+v//RWL/fV6VSyakiAHlijjuAVqulBw8e6Msvv9Rvv/2mb7/9Vn/++ad++eWXvEsDkAM67mD4mCAAgjuAzc1NPX36VG/evNE///yjN2/e6OnTp9rc3My7NAA5YHEyABYngeuHxcngWJwEMIngDqDVaqnZbKrX6+no6Ei9Xk/NZlOtVivv0gDkgE+VBDD+duT6+rr29/dVqVS0ubnJtyaBa4o5bgC4hJjjBoArhOAGgGAIbgAIhuAGgGAIbgAIJvNPlZjZgaRfMz3J9fKNpN/zLgKYgtfmxfqfu5emHcg8uHGxzGxw0keEgDzx2lwcpkoAIBiCGwCCIbjjaeddAHACXpsLwhw3AARDxw0AwRDcC2ZmNTN7YWZbab9gZjvj/YzOuWJmL7J6fkCSzKycdw3XBcG9YO6+J2lLUs3Mau5+KGlH0s8ZnrMr6TCr50dsqZnYNbOt1EiMm4nvZzymbGbPJvZ3JBXMrLGQoq85rsedn4eSnkn6ThOhamZrkgaS6u7eNrMVSRsahXtZ0l7a35JUkzRMj78nacPdD9NjypKGKbSBE7n7npkNJT1NjcQ4iIczHjPU6DU3VkxNyV6WtWKEjjsn6YX/dHKKJHU4g/QXYGhmayl4y+7+3N23J/a7kp5Lepi2dyXV09vVsru3NQp44NzG02zpv9+bWS2N18bTb2msnBoGmdlaOr527Dka6TmOP2cjbe+YWSGvP2skBHeO3H1b0kq6SdJdfei+h2lf+ncXszdl+216zmHq1AsXXS+unymNwoM0vqf0Wh03Gu7evYjmY9F/xogI7vzd02jaQxqF8HiBpyzp5VmfLHU6J85NAudwlumPT24+MB/BvWDpbeXG+C1nmjL5OW1vaLRouSKp5u7b4+3Jt6hpv6xRpz7evptuZY3+4pQ16ngaxx4DTPNCH3e7xfF89xl9cvOB+fgCDgBJ79dYxguSe+4+TP/oP9OoKVjR6B3iPY1CeTxeSNv30kLn9xoF+GTzsTNxfNpzPtRoTeYwNTCYgeAGgGCYKgGAYAhuAAiG4AaAYAhuAAiG4AaAYAhuAAiG4AaAYP4Pks7qXK2ftDUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "#Fix the boxplot\n",
    "bp = ax.boxplot([pe_normal_array, pe_uniform_array], labels=['Normal', 'Uniform'])\n",
    "\n",
    "#We change the axis letter size to see them better\n",
    "plt.rc('xtick', labelsize=20)\n",
    "plt.rc('ytick', labelsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a clear difference between the PE of the uniform point clouds and the normal point clouds. Note that, in this experiment, all barcodes have the same number of bars. Then, we can conclude that connected components from uniform point clouds tend to have a more homogeneous length than the ones from the normal point cloud. It makes sense, doesn't it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we know how to calculate PE using Gudhi. You can use it as a parameter in your scikit-learn pipeline or as a variable in a statistical test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When should I use PE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are plenty of different methods you can use to compare barcodes. Depending on the situation, ones might be more convenient than others. Then, when shall we use PE?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A first remark, is that PE does not keep the moment of birth and death of the bars. Then, if this information is important for your analysis it should not be used. As we will see, the ES function solves this problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fix two barcodes $B_1$ and $B_2$ with persistent entropy $PE_1$ and $PE_2$. In this [paper](https://arxiv.org/abs/1803.08304), a bound for $|PE_1 - PE_2|$ in terms of the bottleneck distance was proven. It can be deduced that, when dealing with  a collection of barcodes $\\{ B_i \\}$, PE tends to be unstable if the number of bars in each barcode is very different from one to others. Actually, in that situation, PE is highly correlated with counting the number of bars. Another cause of instability may be the lack of a lower bound on the total length of the barcodes (in other words, the sum of their bars)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, PE should be used when:\n",
    "- The number of bars is similar between all barcodes    \n",
    "- There exists a lower bound for the sum of all bars in each barcode.\n",
    "\n",
    "In particular, PE is ideal for analyzing 0th-dimensional barcodes coming from alpha/Vietoris-Rips complexes when all point clouds have the same number of points. In this situation, there is a good interpretation of what is measuring persistent entropy. Please, refer to Proposition 6 appearing in this [paper](https://www.mdpi.com/2227-7390/9/15/1723)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Entropy Summary Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is the ES function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we said, PE has the disadvantage that it does not keep information about the birth and death time of the bars. This also happens, for example, when counting the number of bars. In that case, a possible solution is to create a function which count the numbers of bars which are alive at each moment. It is known as the Betti curve and, given a barcode $B = \\{[a_i, b_i)\\}_{i=1\\ldots,n}$, its mathematical definition is the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\beta(t) = \\sum_{i=1\\ldots,n} \\delta_i(t) $$ where $\\delta_i(t)$ is $1$ if $t \\in [a_i, b_i)$ and $0$ otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ES function is constructed based on the same idea. Instead of calculating the entropy of the whole diagram, we only consider the ones which are alive in a concrete moment. Its mathematical expression is the following:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ ES(t) = \\sum_{i=1\\ldots,n} p_i\\log(p_i)\\delta_i(t). $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the ES function in Gudhi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to see some examples, lets generate the 1-dimensional barcodes from the previous point clouds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgms_normal = []\n",
    "for pc in normal_pcs:\n",
    "    rips = gudhi.RipsComplex(points=pc).create_simplex_tree(max_dimension=2)\n",
    "    rips.persistence()\n",
    "    dgms_normal.append(rips.persistence_intervals_in_dimension(1))\n",
    "    \n",
    "dgms_uniform = []\n",
    "for pc in uniform_pcs:\n",
    "    rips = gudhi.RipsComplex(points=pc).create_simplex_tree(max_dimension=2)\n",
    "    rips.persistence()\n",
    "    dgms_uniform.append(rips.persistence_intervals_in_dimension(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:title={'center':'Persistence barcode'}>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAERCAYAAABmRCz9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAP1klEQVR4nO3dTXLjSHrG8ed1+ABsdWvtCtYNGJoTWHUDavoETUX4ACr3CcYsn4DyARzdqoX34pygKe29EB1ea6YaC+9fL5DogiB+gABU4Av+fxGIkvCRTKSgR8kkEmXuLgBAXP/QdwUAAO0Q5AAQHEEOAMER5AAQHEEOAMER5AAQHEEOmdnYzOZm5mb2YGY3aVmY2Z2ZTd7odRdmtniLsr81M7tMbXfXd122MbNR+nlyz/HAGPeRo2BmT5IW7v6ptG4s6UnSe3dfd/x6l5Lk7ssa+96U63WMzGwm6YO7X/Vdl13MzN3d+q4HukOPHFVZ+ZsU3o+SPnb9Qu6+rBPiyZ+6fv038KXvCuA0EeQ4emY277sOwDH7x74rgOOWhlYmkn5K348k/SzpN+W95Ht3X6ZhkrmklaR7ST9K+sXdP5vZjfJe/UjSh7TfKP0rd/+Qyt603zgtk7Q9c/fbPXWZpGMzScUY/AdJT8WxG47PJI3c/fOusms02cjMpunrs3R+5decpPWjarl72nBfXWeSiqGvcWV4bK58eOyLKu+4Wp4rjoW7s7DI3aX8l/1O0lR5MNyk7yeVfUal7x+K79MxT8pDapKWmaRpaf9pUZ6kS+WhoT37TSXdbanvtrpM0/Zx+n6UX+4vjn/YcG6TfWXvaL+p8j8W5XULSbNKOeXz/L2y/6s2rFHXh0pdp8o/65DyPwjl48Yb2uHgc2U5roUeOap+89TTq0o9zczds9LqlaQ/S7pVGiNO2x/TMWeSFunf5baylfcm6+xXty5Zqsu6qJOZycxG6etL5SH/WDr+yt0fa5S9y2Pl+7u0FMddeeVD46JO6dtNbbivrirX1fMe/J2Z/Wf1OHdfm339nLPlueJIEOQ4xFj6erdJcqevb+mlSpB5PtRxLela0tzM1pL+uRIctfc7sC677rIZV7eXAq9O2XWtlfeuC1+KIaJSeWd6OeRR/WOwr66bPmTNlA8nZRu2Vcvu6lzRE4Ich1grf8tde/zUzC7T/sU48EL58MGnhvuNm9Zlw7mMd2xrU3bZSC9D8UGpNy1JRe+40is/tK5nW173N0n/sqd+XZ4resJdK6gtDXdkpTAtJhNd7jhsUtm+bQLQrv3Kvdqxu69r1mVTwBXnspS0Lu9vZhMzmzQ8z0I1cK+VPtQtPugshfiotN/WsmvUdVSp61TSZ3f/r3TcpHxcpew254ojwYQgFL3cqfLAeVR+p8TGyTeVOxykfHy1uGvlo/IgW0i6TWPRs7Rf8fZ/7O6fUqD8rDzAyveov9iv9LoL5b1Z+ea7Vsp1Ke5auUhl/5r2u1E+7jtPY8WjtN9Det3Mv95BsrHsHc1YDE980dcwH0sq/uiUz+NJX4dPRsqHQBbK//C8asNSferU9Ux5D/vThm1Zes371A4f08/o4HPFcSHIASA4hlYAIDiCHACCI8gBIDiCHACCI8gBILheJgT98MMP/u7duz5eGgDCenh4+Ju7n1fX9xLk796902q16uOlASAsM/vfTesZWgGA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiul/9Y4pT8fvOvfVfhle8+/VvfVQDQIXrkABAcQQ4AwRHkABAcQQ4AwRHkABAcQQ4AwRHkABAcQQ4AwTEh6I0x+QbAW6NHDgDB0SPvybFN3eedAxAXPXIACI4gB4DgCHIACI4gB4DgCHIACI4gB4DgCHIACI4gB4DgCHIACI6ZnT1hJiWArtAjB4Dg9vbIzWwsaS5p4e7L0vobSd9L+kXSmaQrd79+q4piu+pzW+jtA6dlZ5Cb2WX6crxll1lalpJ+6rBeAICadgZ50QM3sy8bNmfu/t2b1AoAUFvrMXIzm6ThFwBAD1oFuZlNJa0lTcxsvmffmZmtzGz1/Pzc5mUBACWNg9zdb939s7tn7v5Z0rQ0pr5t/wt3vzg/P2/6sgCAisZBbmaTyqpHSR/aVQcAcKhGQZ5C/K+V1SNJT20rBAA4TKMgd/dHSR8rq8eSfm1dIwDAQfbdRz6RdCnpQtJHMxu7+23avEqTgjJJ75VPCMresK7YgglAwGnbdx/5o/Kx7087tgEAesRDswakOlW/DnrzQHw8NAsAgiPIASA4ghwAgiPIASA4ghwAgiPIASA4ghwAgiPIASA4JgQNCJN7gNNEjxwAgiPIASA4hlYGpsnzVnAcGBpDU/TIASA4ghwAgiPIASA4ghwAgiPIASA4ghwAgiPIASA4ghwAgmNC0MAwqQQ4PfTIASA4euQDdKzT9Hm3ALwNeuQAEBxBDgDBEeQAEBxBDgDBEeQAEBxBDgDBEeQAEBxBDgDBMSFogJh4A5wWeuQAEBxBDgDBMbQyUMf6vJW3xrASThE9cgAIjiAHgOAIcgAIjiAHgOAIcgAIjiAHgOAIcgAIjiAHgOCYEDRQTIwBTgc9cgAIbm+P3MzGkuaSFu6+LK0fSZpJWksaS1q6++Mb1RMdizSFn3cXwG47g9zMLtOX4w2b7yRdu/s67XtvZlfunnVbRQDALjuHVtx9mXrhX8rrU298XIR4spZ0KQDAN9V0jPxCUlZZl0n60KYyAIDDNQ3ykSq9dEl/l3TWqjYAgIO1uWvloNA2s5mZrcxs9fz83OJlAQBlTYM8U94rL/ter3vpf3D3W3e/cPeL8/Pzhi8LAKhqGuQrve6RjyTdt6oNAOBgjYI83WK4SveYFy4kLTcfAQB4K/vuI58ov6XwQtJHMxu7+23afCVpZmZr5b3zn7iHPA4m2QDDsTPI00zNR0mfNmzLNq0HAHxbPGsFAILj6Yc4SKRntNTBEBOGgB45AARHkANAcAQ5AARHkANAcAQ5AARHkANAcAQ5AARHkANAcEwIwkGYQAMcH3rkABAcQY7GhjZdH4iKIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIAeA4AhyAAiOIEdjTNcHjgNBDgDB8dAs9KarKf68M8Cpo0cOAMER5AAQHEEOAMER5AAQHEEOAMER5AAQHEEOAMER5AAQHEEOAMExsxO9YUYm0A165AAQHD1yHJWunr/SJ95p4FujRw4AwRHkABAcQQ4AwRHkABAcQQ4AwRHkABAcQQ4AwRHkABAcE4JwVJhMAxyOHjkABEePHIPR9fR+3h0gitZBbmY3kr6X9IukM0lX7n7dtlwAQD1d9chnaVlK+qmjMgEANXQR5Jm7f9dBOQCABjr7sNPMJmY27qo8AEA9nQS5mU0lrSVNzGy+ZZ+Zma3MbPX8/NzFywIA1EGQu/utu39298zdP0uamtnllv0u3P3i/Py87csCAJLWQW5mk8qqR0kf2pYLAKinVZCnEP9rZfVI0lObcgEA9bW6a8XdH83sY2X1WNKvbcoFmmACD05VF7cfrtKkoEzSe+UTgrIOygUA1NA6yN39Ufm4OACgBzxrBejQIc97YSgIXeHphwAQHEEOAMER5AAQHEEOAMER5AAQHEEOAMER5AAQHEEOAMExIQjoEJN80Ad65AAQHD1yIJhDHgOAt3Us78DokQNAcAQ5AARHkANAcAQ5AARHkANAcAQ5AARHkANAcAQ5AATHhCAgmGOZhILjQY8cAIIjyAEgOIZWAOw05Ge7DGWYih45AARHkANAcAQ5AARHkANAcAQ5AARHkANAcAQ5AARHkANAcEwIArDTUCbNDBk9cgAIjh45gNq6nq5Pb78b9MgBIDiCHACCI8gBIDiCHACCI8gBIDiCHACCI8gBIDiCHACCY0IQgNqYwHOc6JEDQHAEOQAE13poxcxGkmaS1pLGkpbu/ti2XADoSvUZMUMbIupijPxO0rW7ryXJzO7N7Mrdsw7KBgDs0WpoJfXGx0WIJ2tJl23KBQDU13aM/EJSVlmXSfrQslwAQE1th1ZGkr5U1v1d+Vj5C2Y2Uz6WLkn/Z2b/3fK1+/CDpL/1XYlgaLNmaLfD1W+zf5+/bU3ezj9tWtnFGPlZnZ3c/VbSbQev1xszW7n7Rd/1iIQ2a4Z2O9wpt1nboZVMea+87Hu97qUDAN5I2yBf6XWPfCTpvmW5AICaWgV5usVwZWblMfELScs25R6x0ENDPaHNmqHdDneybWbu3q6AlxOCziStmBAEAN9O6yAHAPSLZ60AQHAEOXBEzGxuZsyMxkFObmjFzCbKP5BdK7/DJnP3vR/OmtlU+WcA7yVNJC3c/XNp+1j5c2cWyu/mGSuf4boYwmcGTdqtbps0/Zkcu4Ztdq/tj7hYu/v7oV9rkmRmd5L+Uvd86rT1UK8zSZK7n8yi/IK/q6y7kzTZc9xU+TNliu9Hkp4k3VTKfpL0uyRXfgvmznKjLC3abW+bNC372JcWbbZQ3lEYV5ZZcexQr7V0XgtJ83Rul1219VCvsz/Ope8KfOMLZVG9ONIvzf2e42Yb1k3zv4MvLpRBXBQdttveNmla9rEvLdpsumX9rPT1YK+10jk+HRDke9t6qNdZsZzaGPmflb+tKtv5tMb0NnaR3paVPabt1fVDdHC7HUnZfWp0Xl4ariuY2Y3nj7jAZnXaeqjXmaQT+rAzBfLIXz5yV56em74tkNP+n/T6Ihilf1+sN7ORmU0qk6TCatpulTI2tkkXZR+jLs8rffC5cRx3aNdaE3XaeqjXWdnJBLlePxOmauvDv9z9o7/+jzJ+lPRYWf+jSh+mmNldmjAV2WjP9n0PTdvVJm3LPlajPdsPOa+xb/7Ab4jXWhOjPdvPau4T2ikFeWdKs1mvinXpr/1f3H3p7ln65buX9B/91LJ/tEk7qadYfSdIu+IVgryZO0lX296qlSwlTU+0pySJNmlprvz2wldoV5SdUpBn0h+96U1qPXrXzOaS5v76HtXZht2LMiOPYWZSs3ar0SaNyz5ymdTuvNKxFxsCe8jXWhOZtLet6+wT2skEeeo9Z6qMhxU/3C3jkC+kX6D7DSFe3Nmy7Zfo1dvjKJq2W5026eJncow6Oq9LbbhuhnytNVGnrYd6nZWdTJAnS+X3jpbVeuxuuntgVQ5xMxubWfGfT19Xh1qU3/JU/UA0ooPb7YA2afwzOXJtz+tPm1aewLXWRJ22Hup1luv7RvZvuSh/2/lQWVedATaS9KDSxAzlvaOpXs+2m5f22TT780EDmLjRot32tkmdsiMuTdussu/GySpDvtZK57RxQtCW66xOWw/yOiuWU3zWyqXyv8xr5T/cR3/Zyx5J+h9JH939Nn3/+5bi1u7+vnRsEfZS/kyWub/uOYV0aLuV1u9tk31lR9W0zdK2ufJ7n6+3lD24ay21x8/Kw3qmfNLdUqXhzB3X2d5raKjXmXSCD80CgKE5tTFyABgcghwAgiPIASA4ghwAgiPIASA4ghwAgiPIASA4ghwAgiPIASC4/wcnu2/59HPkcgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#An example of 1 dimensional barcode\n",
    "gd.plot_persistence_barcode(dgms_normal[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate the Betti curves of these barcodes using the [representations.BettiCurve](http://gudhi.gforge.inria.fr/python/latest/representations.html#gudhi.representations.vector_methods.BettiCurve) method. When applied to a list of barcodes, it gives another list with their respective Betti curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample range is the interval where the Betti curve is calculated. It is an optional parameter\n",
    "BC = gd.representations.BettiCurve(sample_range=[0,1.5], resolution = 150)\n",
    "bc_normal = BC.fit_transform(dgms_normal)\n",
    "bc_uniform = BC.fit_transform(dgms_uniform)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, the Betti curve of the upper diagram is the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f78b960b190>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAEACAYAAAB1dVfhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAixUlEQVR4nO3da4wk13Uf8P/p13TvTPcOZx/kUtTscsjQEi1BynKs0BEsyeLQcRIIzgeunA+Jgjy060RJFMCIKAZGjABCbDIfAiVfwpURwAEsy1o6toNIFLzDyJZlObKXI1mSSYo0R6JER8vlPma7Z7cfVd03H6qqX9OP6q5H37r1/wGCuDPdU3emp86cPvfec0UpBSIi0ltm0QMgIqLpGKyJiBKAwZqIKAEYrImIEoDBmogoARisiYgSIBfVFz569Kg6depUVF+eiMg4zz///FWl1LFRn4ssWJ86dQqXLl2K6ssTERlHRF4b9zmWQYiIEoDBmogoARisiYgSwFfNWkROA9gCsAtgTSl1PtJRERHRgKnB2g3UTyilzrj/fl5ELimldiIfHRERAfCXWX8GwJm+fz+ilNqLZjhERDTKxGAtIqsANpRSu26GvaeU2o1lZDSzTkfhxctVNO0OKsUc7j9eju3atYaFV67sAwBOrh3CkZWl2K5NlAbTMutNANdF5DEA2wA2ReRxpdS56IdGs/ryd6/gn/56b237//nF92Pj2Eos1/7kb38bX/j2jwAA73rrKn7vY++N5bpEaTFtNcgqgA0A20qpPaXUNoANN3gfICJnReSSiFx68803Qx4qTXPtVgsAcO59GwCAy9VGbNf+0c063vGWCt7/wDFcvlmP7bpEaTEtWO/BKX3s9X1sF8Cjox6slDqvlNpUSm0eOzZyxyRFqGl3AAB/8/6jAIBq3Y7t2tWGjfW1Q3jgzpVYr0uUFtOC9bj94nshj4NC0LTaAIBjbr241rBiu3atYaFSzKNSzKNutWG1O7FdmygNJgZrN6PeFpGNvg9vAvitKAdF8/Ey62NlJ1hXGzFm1nUblVIelVIeAFCL8dpEaeBn6d5HATwhItcAHAHwONdY68kL1mvLBQBAtR5PZt2yO6hbbZSXcigXc91re+MgouCmBms3u348+qFQUE27jaVcBtmMoLyUQzWmMohXbqmUnDIIgNiuTZQW7A1ikKbVwVLOeUkrpXxsE31euaVSynXLIJxkJAoXg7VBmnYHS/ksAKBcjC+z9sotlWIelZJbBmFmTRQqBmuDeGUQwMusYwrWo8ogMV2bKC0YrA3StPvKIMV8bKtBvJKHk1mzZk0UBQZrgzg1a6cMUinlYltn3ZtgzGG5kEVGuHSPKGwM1gZp2m0s5fsy67jLIMU8RCTWEgxRWjBYG2SwDJJDrWmj01GRX7dat5HNCA4V+ic3mVkThYnB2iBNu4NCtwySh1LAfiv6oFltWCgXcxAR59oxZvVEacFgbZCm1R6YYATiWZVRrVvd63nX5gQjUbgYrA3S6i+DeOudY9icUm3Y3et51+amGKJwMVgbxKlZu2WQGLd9M7Mmih6DtUGcHYy9TTFATGWQxlCw5moQotAxWBtkYAdjN7OOoQxSHyqDFPO41WrDZk9rotAwWBtkoAziBs84NsbUDmTWzrX3m6xbE4WFwdoQSqmBCcaVpXgmGO12B7da7W7ZBehficJgTRQWBmtDeAcPeDXrXDaD5UI28ok+b1u5d+hA/39zkpEoPAzWhugGa7cMAsQz0de/1bz/ugA77xGFicHaEE3bOSzXK4MA8Syh63bcG1UGYWZNFBoGa0M0LS+z7gvWMWxO6WXWg5tiANasicLEYG2IXs26rwwSS2bdO3ige132tCYKHYO1IUaWQUoxBOvGwWC9UshBhDVrojAxWBuiN8HYX7PORX4IgPf1+8sgme7p6iyDEIUlN+0BIvIJAEcA/BaANQBnlFLnoh4YzaZXsz64GkQp1W1fGrZq3UJGgOXC4K9SHFk9UZr4zazPAngOwDkAj0c3HJpXtwyS772k5WIOHQXcarUju261YWNlKYdMZvCPQbmY5wQjUYimZtYA9pRSd0Q+EgpkdBmkt97Z29EYtmrdGqhX966dY2ZNFCLfd7CInIYTuHcjHA/NaWSwdoPoH3z3TTxw5wpOr99xIAP2vHS5iivV5sRrvHt9tfsH4PLNBl5+o4bvX7s1sCGm/9rfvVzDV15+EyePHMLJI8tzfV9E5PAVrEXkMQDbALZE5JxSamQpRETOwimZYH19PbRB0nRNy1sN0qtZ33W4CAD4d7/zbQDA584+jIc3jhx47q2mjQ/916/Cak8+r/EfPLyOT/29dwIA/tVv7uDPvn8DAPDBtx0/8NgTh4u4+MIb+Mh//1OcOFzEnzzxyBzfFRF5pgZrpdT5vn8+IyJPishFpdT2mMeeB4DNzc3oT2qlrlGZ9V9/6yqe/fhP4eU3avj4576JK7XRmfON2y1YbYWP/fR9IwMvAPzi5/98IPN+o9rE+x44ho8/cj/uP1Y+8Pgn/vbb8XPvvhuf/foP8bvf/KtIJzmJ0sDPapDTSqmdvg/tAHgUTqZNmmiN6A0iInj7iQqOLBcAjG+X6i2/e+dbDuOhk2sjH3O8XBxYBlhrWDh15NDYx5cKWTx0cg2Xvn8D7Y5C3WrjUCGaujlRGkxcDeLWqZ8b+vAqgFejGhDNZ7jrXr9eY6XRqzO8zSvlEbXn3tfoTRgqpZxzFyc83u+1icificHazaiH69MbAD4f2YhoLt7SvUL24Eu6lMsgn5WxqzOq3Y0t44NvuW/r+u1WG+2OGmiLOv55bJdKFAY/70svuRtj9gDcB2dTzF6Ug6LZNe0OCtnMyNUeIuL0CRmz/bvX32P8r0Ol2GsKNWqL+fjnsV0qURj8TDDuwKlTk8aaVmdgcnGYs6NwTBlkRE/qUc+vNSx0OqrXFnWWMggza6JA2BvEEE27PbJe7XEy43GZ9cHTXg4+P+/uhrT7Muvpb8y8niGsWRMFw2BtiP7DckeZ1Kuj1rCwXMgiN6Le3Xu+V3u2e2UTZtZEsWGwNkTTnlIGmVSzbozeMj78fMCpPXc77fmoWXcnGFmzJgqEwdoQTauNwsSa9fh2qdW6PXVlR/+5iqNOhxlnKZdFMZ+JvFUrkekYrA3RtDsDp8QMm3RqTLVhTS1p9M5VtH2ty/Z7bSLyh8HaEE27PbEMUi7m0LA63fXY/fyUQfrLGdWGjWI+MzGTH34uJxiJgmGwNsTUmrUbjEeVI6p1e2pJo/d8y2mL6jOr9p7LzJooGAZrQzjrrCeXQYDRE30zZdYN29fjh6/NCUaiYBisDTF1nXXf0rt+SinUfPT5yGczOFTIOmUQH5n44LXHb8ghIn8YrA3RtDtYmrROekxmPUufD2+isDZzZj1+Qw4R+cNgbQhnNcj0mvVw7XimPh8lZ6LQb8e9/mtXG87BvUQ0HwZrQzStts+a9WA5YqY+H25m7Zy7OEMZpJiH1VZouCewE9HsGKwNMX01iBNchw8gmKnPRymPm+6mmNky69HXJiL/GKwNoJRCqz05WJfyWWQzB3taz9Lno1zM4UqtCautfG+IcZ7H/iBEQTFYG8BqKyiFiTsYnZ7WBzenzNqb+k33HMfZyiDOY29yYwzR3BisDeDtSpyUWQOjN6d0mzL5WQ3SF6BnnWAEmFkTBcFgbYBRJ5uPMmpzyix9PvoD9KybYvqvRUSzY7A2QHPEyeajOIfeDpdB/Pf56A/Qs22KGb0hh4j8Y7A2QNNyyyAT1lkD4zNrvyUNZtZEi8NgbYCZyiAjlu75Dbzz1qyL+SwKuQxr1kQBMFgbYJYyyHDXvVn6fPQHaD/b04efywMIiObHYG2AbhnER2Z9u9WG1e7tJJylz4f3uKVcBsUJywRHP5f9QYiCmClYi8iWiDwW1WBoPt3MekrN2suG+zPcWfp8eM+fZUNM77nsvEcUhO9gLSKrAJ4GsBbZaGgu/ssgByf6qnXLd0nDe9wsG2K612bnPaJAZsmsPwxgO6qB0Px8b4oZ2vatlJppgtE7/HaWycXutXlaDFEgvlIkEdmCE6gfinY4NI+mNVtmvf3iFbx+ow6r3YHVVrPtRizmZ1q21/+8q7UmvvjtH419zNpyAQ9vHJn5axOlwdRg7ZY/VpVSuyIy7bFnAZwFgPX19TDGRz74rVmfOFwEAPyX514Z+Pjdq0Xf1zp1ZBmnjhyacYTAPXeUUG3Y+Be/sTPxcV/75Adx92pp5q9PZDo/mfWWUuoZP19MKXUewHkA2NzcZKf5mHhlkMKEk2IA4K1rh/C1T35wYIIxnxXce3TZ97V+/Z+8B9nM5D/ao/zC++/Dzzx4Jzpjfiv+7+41/PL/+gtc228xWBONMDFYi8hpAJNTIVo4v5k1gMCBsFSYbcmeJ5sR/LU7y2M/v3e7BYDNnojGmZZZrwHY7Ct/bAFYExEviyYNeDXraZm1zsrckk400cRgrZQaWP0hIo8CuMhArZdWu41cRpBLcLDuNXtisCYaxfeCWXfycAvAqohc91vHpug1rcmnxCRBbw04N84QjeI7WPdPHpJenJPN56sl62KlkIMIM2uicZKdjhEAZzVI0jPrTEZQXjrYaIqIHMm+wwnA9JPNk6JSOthvm4gcyb/Dya1ZJ7sMAozut01EDgZrAzTttq811rpz2qiyDEI0SvLvcDKmDFJmZk00VvLvcHKDtSFlENasiUZisDaACatBgNGnrxORI/l3ODkTjCbUrIt57DdttMd1eyJKseTf4WROGcTdxbjP7JroAAZrAxhTBimyPwjROMm/w8mY1SBeZn2Tk4xEByT/Dic0rQ4KJgTroTMiiagn+Xd4yiml3DJI8mvW3unp3BhDdBCDdcLZHYWOmn6yeRIcLjGzJhon+Xd4ys1ypJfuvDIIO+8RHZT8OzzlWl6wNqAMstItgzCzJhrGYJ1w3snmJpRBsm5Pa5ZBiA5K/h2ect5huSaUQQCvpzXLIETDzLjDU6xpUBkEcFaEMLMmOojBOuFMKoMAPC2GaBwz7vAUMy2zrhTzXA1CNAKDdcIZV7NmGYRopNy0B4jIaQBrAFYBbACAUuqpaIdFfrEMQpQOU4M1gOcAPKKU2gYAEVEisq2U2ol2aOSHeWWQHGpNG52OQiYjix4OkTb8pGOPeIFZRFbdj+1FNSCajYmZtVLAfot1a6J+UzProQz6wwCeUUrtRjckmoV5NWtny/nn/vQHuOeOQ/jZH79raob9/GvX8eCJwygV9Hh30bI7+NJfXEaj1fb9nGxG8OiP39n9/ud1+WYD+00L9x8vB/o6pB8/ZRCIyAaALQCPKqXOTHjcWQBnAWB9fT2UAdJkppVB1o8cAgD8xy++BAD4vY+9F+966+rYx9+8beHMf/sT/Iefewf+4cMn4xjiVH/8l1fxr3/zGzM/75duvx3/7Kc2Al37V599ES9druFL/+Z9gb4O6cdXsHYz6fMiAhG5MC5gK6XOAzgPAJubmzxILwamlUEe3jiCS7+0hW/8YA8f/R+XcP12a+Ljb9xuoaOAq7VmTCOc7totZ8zP/MJP4sRqyddzPvCfvozrtyZ/r35c3W/h6r4+PwsKj69g7VFKnReRJ0XkE1wRooduGcSQYA0AR1eWcO/RZQDTmzp5y/x0Wu7njfn+4ytYPVTw9ZxyMR/K91BtWKjWbSilIMIJWpNMvMNF5LSIvDr04V0A90U3JJpF0+4gI0Aua06wBoBKyTuPcfJEo9dHRKd+Il7QXVnynwtVirlQvodq3UKr3emWx8gc0+7wPQDbQx/bAHAxktHQzEw5JWZY94ivRGbWNlaWcjP9Aa2UwsqsvT9e+vw8KBwT//QrpXZF5II7cQgADwF4XCn1TPRDIz+adseYlSD9ivksCrnM1ADmBSWdglO1YXVPaverUgy+GUgp1ft5NCwcrxQDfT3Si5+le8OZNWmkaZlxsvkoTgCbUgbpZtYalUHqVvekdr8qpRwuVxuBrlu32rA7zrz+TY3KQhQOM+/yFGm1O0aWQQAngE3PrPV72+9k1jMG6xAy6/4/bDqVhSgcDNYJ59SszXwZ/XTgq7lBqaZRcKo17O4EqV+VUvBug/0/A3YuNI+Zd3mKNC0za9aAv6ZOXvnD6yeig/ky6xzqVrt7pua81+3+t0bvNCgcZt7lKdK0zS2D+Dk1xgtKOvUTqdZtlGecYCx3T3afP8iyDGI2BuuEM70M4neCEdAjm+x0FGqN+SYYgWATpYM/Cz3+cFF4zLzLU8TJrM18Gf1OMGbdRk86BKhbLRsdhbkmGIFgf3C852YzwszaQGbe5SniLN0zswxSKebRsjtoWOO711UbFu5y1xPrEKC8zHieCUbn+QGCtXvtuypFLd5lULgYrBOuabeNnmAEJgewat3CW+4odf970bwxzJ9ZByiD1C0s5TI4ulLQat05hcPMuzxFjC6DuJN04wKY3e7gVquNe7xgrUGA6gbruWvWQTJrp1bOo9HMZOZdniImrwbxAt64FRL7TSc433PHoYmPi5O3vnnezDrQapCGjUox565PX/zPgsLFYJ1wTcvs1SDA+IzZy7jvWS0N/HuRvMx41pr1oULWmRgMWAZxMuucFu8yKFxm3uUpYmojJ6C/DDI6S/QC4x3LBSwXsnpMMLpjLc+YWYuIr3XlE6/dsFEp5kPZuk76MfMuTwm73YHdUcaXQcYFsN5kXk6bOq2X0c66KQYI3h+kVu/VrJtTVtFQ8jBYJ1irbd4pMf2mrZDwgnjZyyY1yawPFbLIz3EYRNDyRbVhoVzMdd+RsD+IWcy8y1PCO9KrYGiwLuYzyGfHb/DwgnillHMCnSY163lPKA+SWTu9rJ0ySLkYfM026cfMuzwlTDvZfJiITAxgvck8nTLr2TvueYJ8D027g1a70/3D5Yxl8T8PCg+DdYKZdrL5KM5xV+NWg1gQAVYKudCOxQoqUGYd4N1B/2acaatoKJnMvctTwGunaepqEMA7SHZcZm2jvJRDJiOhHTgblNPLOkAZZM4/OAPvMqasT6dkMvcuTwHTyyCA15R/fBnEC0xldyOIUovtaT3P+YueSimP26027PbsPa27PUncTTGAHuvOKTwM1gmWijJIcVIZxO4Gpkoph44CbrUWu1ytWrdmXmPtCbKKo3+bexhb10k/5t7lKeCtBjE5WJcnlkGsbmAKo8VoUEopZ2PKnBOMQVZxVPu2uZfyWeQywglGw5h7l6dAtwySN7sMMmlTTLmbWS9+udrtVhvtjgowwTh/+aJ/g5CIaDPhSuGZmgKIyGkAW+4/fwLA00qp7UhHRb6kowySQ8PquCfiDP5RqjX6yiAa1Gn7J/nm0d1eP1dmPXjtsiYTrhQeP+/XtpRSTwGAiKwC+J6IPKKU2ol0ZDRVb4LR4GDdXdlgY2llMFg7jYvcMogGa4u7m3QCZ9ZzBOu6jUI20/1d0GXdOYVn4l3uZtVPeP9WSu0BuIRepk0L1K1Zm1wGGVOLbncUas0RmfUCA9S8Hfc8QUo5Xv1eRLpjYM3aLBN/q5RSOyJyZujDGwD2IhtRir3yRg1feeWq78c//9p1AKZn1s6v6Ge//gOccFuhAr0SkBfgvP+/+MIbuHHbCVIC4GffcRfu7ntelLwlhvNvN3e+1+devIL95myrWr71+t7AdSvFPK5Um3ONg/Q0NQXor0+LyAaANQCfH/VYETkL4CwArK+vhzTE9HjySy9h+8UrMz1nbbkwV4e3pDh5ZBn5rODXvvq9A58TATaOLQNwAt3x8hKe/c5lPPudy93HvH6jjn//oQdjGWuvV8l8wXq5kMOJw0X8/gtv4PdfeGPm5//Mg3d2/5tlEPPMepc/DeARtxxygFLqPIDzALC5ubnY3QkJdOO2hb9x7xrOf2TT93OK+YzRm2LuO7aCb/3y3+p2GOyXywiWl5xf4Vw2g68+/kHU+9qC/p1P/xH2brdiG2uvC+B8fzwzGcEf/tufHvgeZlFe6l1Xl8ZWFB7fv1Ui8gkAT3JiMTrVuoX7j6/g8JyZmalKhSxKmP4HqZDLDHQgPBzz8rXewQPzv9MZ/h7mVSnmUbfaaNkdY7sypo2vV1FEHgOw7ZVE3HIIhSxIEyA6KO7sstqwtXmnU+7uhmQpxBRTg7WIbAHY8zJqd/ne6YjHlUpB2mvSQXHXbat1ff7Y9laWsBRiiomRwc2gL7r/3f+phyIcUypZ7Q7qVlubm90EcR/11d9YatF02H5P4Zq2dG8XzgooipjXvEeXm90Ek5pARcFpLKXHOyMdtt9TuDjzoIkwJqdoUKWUw37TRrsTz8IkrTLr7o5OlkFMwWCtiWrADRV0kNfkaT+m7Lq/V8mieePgBKM5GKw1EXRDBR0UpDHSPJwugHq8M2IZxDwM1poI2leCDvIC1s0YJtmcXtb6lEGWC1lkhGUQkzBYa6L/wFMKR5zNnRpWB1Z7/l7WYRMRlLnl3CgM1prgapDwxTnJpuM7I3beMwuDtSaqDQsZcd6+UjjizKx1fGcU99JFihaDtSa8I6qGNh9RAEGa+c8q6CkxUagU490URNFisNZEkINWabSVpflPC59V75QYfV7DSinHmrVBGKw1oVNfCVNkM4LyUjwBS9fMOo4/VBQPBmtNsONeNJz+IHFMMAY7fzEKcfdGoWgxWGuCHfeiUS7GlFlr2C6gUszjVqsNe8TBDZQ8DNaaqDGzjkRc2WW1YaGQy6Co0eHFvZ7WLIWYgMFaE9WG3e1lQeGJa/ma03FPr9ePW87NwmCtAbvdwX6TZZAoxLUxxNlqrtfr1+2Nwi3nRmCw1sB+U7/JKVPEdVqMjqt5mFmbhcFaA+y4F51K0elp3Ym4p7WzTl6v14+nxZiFwVoDvV7Wer2NNkGllIdSwH4r2lKAM0Gs1+vnlWU4wWgGBmsNdPtKaJaZmSCu7NJZeqnX68cyiFkYrDWg44YKU8TVea/a0OfgAc9KIQcRlkFMwWCtAS/z0e1mN0EcnfcaVhstu6PdH9tMRrCylGPnPUMwWGuAZZDoxNF5T8e+IB523jPH1GAtIhsickFEtuIYUBpVGzZEgPISM+uw9TLr6LJLHTvueSolnhZjiom/XX0BeiOGsaRWtW5hZSmHTIa9rMNW7m4MSWtmneOmGENMDNZKqW0AEJHr8QwnndhxLzpx9MfQ8ZQYT6WUxw+v3170MCgE+r1vS5E/+/51fO0vr+Fbr9/UMiszQS6bwXIhO7EU8M0f7uEPv/vm3Nd4+UoNAHBYs+3mgPMH5P/t1fHp7VfGPubOyhL+/nvWYxwVzSPU3y4ROQvgLACsr/PFn+ZT//sF/PnrNwEAZx66Z8GjMde0zntPfeklfO3Va4GucWS5gLsOlwJ9jSi84y0V/PbO6/jP2y9PfNwHfuw47jpcjGlUNI9Qg7VS6jyA8wCwubkZ7f5eA+zVLXzoXXfj0z//bvDoxehM6w+yd9vCI287js98ZHPua4hAy/Mz//F778U/+slTYz//xe/8CP/ys9/AzbrFYK05/d63pUi1buFwiROLUXM6742vWVcbFt52V9nY12HS93WYuxwTg+usF0Qp5TT/0XBSyjTTMutq3UrtnAGbPSUHg/WC3G610e6o1AaJOE062qvTUag1bS3XSMeB/UOSY9o669MAtgBsAnhcRDbcujQFVGM/kNhMOjR3v2VDKT3XSMeBBxQkx7R11jsAdgA8Fc9w0qO3kSKdGV2cKsU8ag0LSqkDk4DeH8209mXxjpKrMbPWHssgC6LzRgrTVEo5dBRwq9U+8Lm0vw6FXAalfJbNnhKAwXpB2GkvPpMm0dhEy63pc4JRewzWC8KjvOIzaRKNvcTZ7CkpGKwXpHeUV3qDRFx6mfXBt/q9zDq973DY7CkZGKwXJO0TW3Ga1HmPfzSZWScFg/WCVOsWlnIZFPPZRQ/FeBPLIHX+0eQBBcnAYL0g1UZ6d83FrTIlsz5UyCKXTe+tUCnx6K8kSO9v6IJV6+ndNRe38oTTYmrsJT6wDp30xWC9IMys4+OtJR618aNat1M9uQg4ZSKrrdCwOoseCk3AYL0g1brVzfgoeuM67/GUnr4JWE4yao3BekGcjnvpzujiNK7zHt/hsPNeUjBYL0iNQSJW4zrvce6AnfeSgsF6AZRSbpBgsI7LuM57zKzZeS8pGKwXoGl30Gp3Uj+xFadRZRDnj6aV6jXWADPrpGCwXoC0d3pbBGeCcTAY3Wq10VF8HVizTgYG6wXo9bJOd5CIk5NZ2wNridlxz9FbDcIyiM4YrBfgptdxL+Vvv+NUKeXR7ijUrV5Pa57W4yjms1jKZVgG0RyD9QL0elmnO0jEaVTnPZ7W01Mujj/6jPTAYL0AXkZ3mEEiNl5A7s8eOXfQ4/QHYWatMwbrBWCQiF95xCQaT+vpYec9/TFYLwAnGONXGbGlmqf19Dg9rVkG0RmD9QJU6zYK2QyWcvzxx6W7lri/Zl1nZu2pFHOoMbPW2tTfUhFZBXAWwC6ADQDbSqmdiMdlNGfXXA4isuihpEZ3grExWAYp5jNYyvEACJ4Woz8/KcUFAOeUUrsAICIXReSMUmov0pEZrFpnp7e4jTrai1v+eypcDaK9ie/D3ax6wwvUrl0AW1EOynS1ho0y66Sx8tYS1/rqsrUm+4J4KqUcWu0OGn3r0Ekv0zLrTQB7Qx/bA/AogGfCHky7o/Crz74Y9pfVzkuXq3jgzvKih5E6lVIeX3nlKjpfeAEA8K3Xb+J4eWnBo9KD9w7jU194ASWeCxrI2vIS/vkH7gv9604L1qsArg997Bqc2vUBInIWTn0b6+vrMw9GKYXf+PoPZn5eEm2eXFv0EFLnPfeu4csvXcFr1251P/Z333ligSPSx4N3V7B6KI//ufNXix5K4p06shxJsJZJ566JyGMAnlBKPdT3sU8A+Aml1JlJX3hzc1NdunQptIESEZlORJ5XSm2O+ty0tWN7cLLrfkdwMNsmIqIITQvWlwAMv19fBXAxktEQEdFIE4O1uzzvkoj016g3AWxHOSgiIhrkZ531GQBnRWQXTpb9Ua6xJiKK19Rg7Qbmp6IfChERjcPmFERECcBgTUSUAAzWREQJMHFTTKAvLPImgNfmfPpRAFdDHA6Fg6+Lnvi66Gme1+WkUurYqE9EFqyDEJFL43bx0OLwddETXxc9hf26sAxCRJQADNZERAmga7A+v+gB0Eh8XfTE10VPob4uWtasiYhokK6ZNRER9WGwJiJKAD+NnEIjIqfhdO3bhdNqdU8pNbWD37zPI3/m+fm6nRgvAHgaTivdDTjHvT2tlNqJdMApIiIXAPyK358p75V4zPK6hHavKKVi+Z87wAtDH7sA4HQUz+P/YnldXgVwA4CC0+Ocr0l4r8nTAJ50f75bUb6W/F8sr0vgeyXOb/Tp4W8OwGkAF6N4Hv8X+euywSAQy+vz6gxBgfeKnq9LKPdKnDXrD8N5a9ZvF8BWRM8jf/jzNQdfS4PFEqzdms2qUmrgF0m5hxi4dbbQnkf+hPHzFZFVETk9dJoQxYz3iv6C3itxZdarUz4/fM5j0OeRP6tTPj/t5/vz6JvMEpELIjLta1I0Vqd8nvfKYgW+V7h0j+biZnC/opTaVkrtKWdW+yKAzyx4aERaCeteYbCmuamDZ3FuA3iM2TXRoDDulbiC9R7g1GzGfP56yM8jf/aA+X6+InJ2wuNZv47fHsB7RUdh3SuxBGv3bcAehupm3i+WGrMwfN7nkT/z/nzdCZKnJ0yUDK9IoIjxXtFTmPdKnGWQbThrPvttuh+P4nnkz8w/XzcwnBteeQBn6djOiLd8FA/eK5oJ816JM1g/DuCJoY+dcz8OoLu05XkReWyW51Eg874u1/uzBTeDOwfgoxGOlVy8V/QU5b0Sa4tUEdmC85d/F06tZkf19S1wv4nvAXhcKXXe7/MomACvy2Po1dzuA/DkiAyCZuT+vJ+AsxzvLIAdONnxRe914b0Sv4CvS+B7hf2siYgSgEv3iIgSgMGaiCgBGKyJiBKAwZqIKAEYrImIEoDBmogoARisiYgSgMGaiCgBGKyJiBKAwZqIKAH+P5wrKafQs8F8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xseq = np.array(range(0,150))/100\n",
    "plt.plot(xseq, bc_normal[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We proceed similarly to calculate the ES function. We use the [representations.Entropy](http://gudhi.gforge.inria.fr/python/latest/representations.html#gudhi.representations.vector_methods.Entropy) method. This time, we include the parameter mode='vector' to specify we do not want to calculate PE but ES. Using a list of barcodes as an input, we obtain a list of curves, each of which being the corresponding ES function. You can use the parameter sample_range to specify the domain of the function and resolution to specify the number of samples you use to calculate it. The normalized parameter indicates if you want to normalize the ES function by its L1 norm (it is done by default). It is used when you want your function to be independent of the coordinates used to generate the persistence barcode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "ES = gd.representations.Entropy(mode='vector', sample_range=[0,1.5], resolution = 151, normalized = False)\n",
    "es_normal = ES.fit_transform(dgms_normal)\n",
    "es_uniform = ES.fit_transform(dgms_uniform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f78b95a58e0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjC0lEQVR4nO3dfZBb13ke8OcAFxdYALveXe6KliVS5NKSkjh2HXrtVLZTpwk1TpvONEnFpl/xTBKLjDOt22kmZt1m2mna1JXb/tFMpx0pk86kjepUVJo4TtK6pJ1p0lhWvaSs2Ioi2VzqW6SWXGIJ4C4+LnD6x70HwAK4HwAugPvx/GY4EoGL3cvF4sGL95x7jpBSgoiI4iM17xMgIqJgMdiJiGKGwU5EFDMMdiKimGGwExHFDIOdiChmtHmfAACsra3JY8eOzfs0iIgi5dKlSzeklOv9t4ci2I8dO4atra15nwYRUaQIIV4edjtbMUREMcNgJyKKGQY7EVHMMNiJiGKGwU5EFDMMdiKimAnFdEdKlld3DexU6tBSAt915xK0NOsLoiAx2GmmjIaJH/x3/weNVhsA8It/9V342APH5ntSRDHDUolmqmQ00Wi18VMfOg4hgJ1yfd6nRBQ7DHaaqWrdBAC89+gyCrqGar015zMiih+2YhLuqSs38c9+55tomG20pESrJbG+lMP5sw9A14J/3682rCAvZtPI62kYDTPw70GUdKzYE+6p7Zt48XoF77l7Ge+/ZxV3Li/g2VdL2K02pvL9VMWe1zUUsxoqdQY7UdBYsSdcpWaimNXwy3/zewAAv3npNVx6+RYaZnsq308FezGrIZ9Nw2iwFUMUNFbsCVepN1HMdt/fVful0ZpO4Kogz+tp5HWtE/REFBwGe8JV6y0UsunO31Ww15rTqdgrPRV7QWfFTjQNDPaEK9dNFHOZzt+7Fft0gl0NluazGgpZDVUOnhIFjsGecJVaE4s9rZisCvYp9dgr9vTGfCZtT3dksBMFjcGecJW6eaDHPu1gN+om8noaqZSwBk85j50ocAz2hLN67D2Dp2mr316f1qyYhom8bn2/gm61YqSUU/leREnFYE+4cq2JxdyQWTFTm+7YQtEerM1n02jL6b2JECUVgz3BpJTOrZipTXfsVuzq+7LPThQsBnuC7TdbaEugOMOKvfeNRAU814shChaDPcHUnPLCkAuUptUeMRot5O1WTEG3/sspj0TBYrAnWKVmBepspzuanTeSvP1fLgRGFCwGe4L1XgWqTL1ir7c6lboaRGUrhihYDPYEUxX7gR57etqzYrqDp90eOyt2oiAx2BNsWMUuhICeTk2lYpdSotroDp4WVLBzvRiiQHkGuxBiQwhxXghxyuO4TwkhHhFCnBRCnBJCPBrcadI0DAt2wOqzT6NirzXbaEt0Bk/Vf9ljJwqW63rsPWG+4fPrnbH/XATw8ATnRTPQCfbcwV8DXUtNZR67mv0yULGzx04UKNdgl1JeBAAhxK6Pr1WSUq4EclY0E+Xa8Ipd11KoT2HZXrUujOqt5zIppAQrdqKgBd5jt1sxfit8mqNq3YSWEp0pjopVsQcf7N3Wj9WCEUKgoHN7PKKgBRrsQoiHAGwDOCmEeMTj2DNCiC0hxNbOzk6Qp0E+VeomijkNQogDt0+rx95Zi13vfkLgCo9EwQss2KWUj0kpn5RSlqSUTwJ4yG3A1T5+U0q5ub6+HtRp0AjUfqf99CkFe/dK1+6OTWqFRyIKTmDBLoQ42XfTZQAPBvX1KXjlukOwT2m6o9oGr3cJA25oTRS8QILdDvUv9d28DOBKEF+fpqPqEOxZLe2rYh91HfVOxd7TimGPnSh4Ywe7Pb99AwCklJcBnOs7ZAPAExOcG02Z6rH307UU6h6Dp99+q4KT/+ICnn/ztu/vZwxZdKyQ1TgrhihgrsFuz3D5FIBNAOeEEGd67j6Lg2G+ZV+kdMYeOD0tpSwFfsYUmEl67FdvVHHLaOLxp1/2/f3UFaZ5vdtjz+scPCUKmtc89suweuWfHXLfOYdjKSIce+xaCnXTPWxVlf35r7+BX/jh70Iuk3Y9Hhg+vZKDp0TB41oxCebcY/eu2Pft6rtcM/HF5675+n5Go4W8nj4wvbKQ1VixEwWMwZ5QrbaE0WgN7bH7CXbVVlkr6nhi61Vf37N3LXalkE1zQ2uigLm2Yii+nBYAA/xNd9y32yd/6wNH8ctf/jb+xmNPQcCqxHOZFP7lj74bdy0vHHiM0RgM9ryuoS2tBcIWdO92DhF5Y8WeUK7B7qNiNxotaCmBn3jgGD5y3zrabetTQLVh4g9e2MGll28N+Z7dTTYUdbES++xEwWHFnlBVh5UdAXseu8d0R6PRwoKexvpiFr/2Ux/o3P5GaR8f/Ndf7kxtPPCYIa0YtbyAUW8BxZH/GUQ0BCv2hHJa2RGwKvZWW8J0Cfd9eyC0nwruYZtnVHp2T1KKrNiJAsdgTyivVgwA16rdaLYGQhrozlEftt2d0Wh1grx7PLfHIwoagz2hhu13qqh55m59dqNuDq3YM+kUdC01tAI3GibyQ2bFANwejyhIDPaEqvqp2N2C3aEVAwAFh6tJK0PmzXd77KzYiYLCYE+osh2ki9nMwH162vq1cJvyaDRbWBjSigGsPnt/xd5qS9Sa7YE3g6JLT56IxsNgTyjViilkB6tuVbG7Bft+w0TeYRmBgj54NakK+oLeX7E79+SJaDwM9oTa22+ioKehpQd/BbKaFbbjtmLy9tWkvapDVnbs/ft/e/oV/Ozjl/DsqyXf/wYiGo7BnlB7+028bWGwDQP0DJ56THd0ulK0oGsDFbiaXrnYN1ib1VL4y+9+OwDgf33zGr7w7Bv+/gFE5IgXKCXU3n4TSw7B3mnFNJ373tUhywMoeT2NG5X6gdvKtSaAwWAXQuA//u33AQA+8EsXO28ARDQ+VuwJtbffwHLePdidKva2PRC64NBjLw4ZPL1tB7bTm4m6r1xvep47EbljsCeUr1aMQ499vzm4YUavfHZwuuPtfSuwl4bMm1cWcxpu77NiJ5oUgz2h3ILdax67MWQnpF7D9jHt9tidK/bFXKbTsiGi8THYE6pkNLGc14fe5zWPXW2y4TSPPa9rqJvtA2vNOA2e9lrMaeyxEwWAwZ5AtWYLdbPt3IrJuE93NJpW+DpW7PbceKNn8LVcayKdEo59ecBq09xmsBNNjMGeQKrf7diKURW7w+Bpta4qdqdg71mK11aumVjKaQe2xevHVgxRMBjsCVTyCnavwVO7FdN/FanSuZq0Z2ZMudZ07a8DVsVeN9ueG2kTkTsGewLteQR7trOkwPCANRoerRh9sGK/XTNd++tAd2CVfXaiyTDYE2jPsILdcR572t90R6dWTN7usffOjLEqdq9g1+xjGexEk2CwJ5BXKyaVEsikxdjTHdWKjcaBVozp2YrpVuzssxNNgsGeQF6tGMCq2j2DPeM83RE4uBSvNXjqFeys2ImCwGBPoD2jASHcLxbStZTjPHa1KYbzrJj0geMA4LaPVswSK3aiQDDYE2hvv4nFrIZ0ynnqYVZLu8xjb0FLic7smX79FXu7LVGpm67LCQDdip1z2Ykmw2BPoL1956tOFV1LOS4Ctu+yFjtgbY0HdNdgrzRMSOn+CQHoVuxqnj0RjYfBnkAll3ViFKsV4zzdMe8whx0AtHQK2Z4Nrf0sJwB0N9Zmj51oMgz2BHJbAEzJau6Dp24VO2BdfarmsXfXYnf/numUQEFPM9iJJsRgT6C9/Sbe5jCHXXEbPHXbPUnJ6+mBin1pwXtfl6UFLitANCkGewLtGT5aMS7THasN07ti10ev2K1juMIj0aQY7AkjpfTVivGu2N2r794Nrf322K1juIsS0aQY7AljNFow2xLLnj12l+mOjVZn5ouTYra7ofXtkYKduygRTYrBnjBeywkoWZfpjobPHru6QrW7LZ6fVgx77ESTYrAnjNcCYIruMitmv+ljVoyuHWjF6PYUSC/ssRNNjsGeMGqdmCUfg6fjzmMHDm5orVZ2dNtkQ1nKZRjsRBNisCfM3n4DgI9WTGZ4xd5qS9Sabdct7gBrHrtatrfsYy12ZTGnodFqo9bkZhtE42KwJ4yq2D2XFHCY7qjWYvfTilEbWvvZPUlZ6qwXwz470bgY7AnjZ8lewHm6o9fuSYq632i2RqzYuYsS0aQY7Amzt9/sXLrvRtdSMNsS7bY8cPt+Z5MN96Du3dD6dq3pa0YMwDXZiYLAYE+YSs1EMes9kJnVrODvn/LotXuS0ruh9SgVuxrU5ZRHovEx2BOm2mh1tq5zo3c2tB4e7F7z2Is9FbufbfEUVuxEk/NXRlFsVOve67wA3WD/lT/c7mxODQCv3DQAeLdi1P3lWhOV+jg9dlbsRONisCdMtdFC3kfFvrFWQDol8B/+4NsD9y1k0rh7ZcH18Wp7vGdeLQHwt5xA73FcVoBofAz2hKnWTRSz3hX7h965hud/8YfQlnLgPi0loKXdu3gr9nTKf/PFFwAA64tZX+dX1DUI0Z29Q0SjY7AnTLVuYrWQ93Ws056mfhxZzeP8zzyAW9UGdC2FD55Y8/W4VErgnetF/Or/vYqN9QJ+7OTdY58DUVJ5BrsQYgPAIwAelVJedDluGcAZANsANgBclFJeDug8KSCGz8HTILz/2OpYj/v1j38vPvm5Z/APn3gW//wLf4phE3jO/oUT+MT3n5jwDIniyfUVLoQ4Zf/vho+vdR7AWSnltv3YC0KI01LK0mSnSEHyO3g6T4eXcnj849+L//rVl/HSjerA/b/3jTfxtZd28Qkw2ImGcQ12VaELIXbdjrOr9Q0V6rZtAKcAPDnhOVKAqg1zZhX7JLR0Cj/5oeND73v+Wrmz1jsRDQpqHvsmgFLfbSUADwb09SkAZquNWrPtOVUx7Ao9+6kS0aCggn0ZQH9VfxOAY5NVCHFGCLElhNja2dkJ6DTIjWEv4FXwMSsmzArZ7n6qRDQoyCtPRxopk1I+JqXclFJurq+vB3ga5ES1LwoRaMW4KfYsCUxEg4IK9hKsqr3XIQxW8TRH1bqq2KMd7HldY4+dyEVQwb6FwYp9GcCFgL4+BaBTsYd8VoyXYjYNo9kaWHmSiCxjB7sQYsOe4w57SuOW+rttE4DjvHeavWpnLfWIV+xZDVJ2N/0gooO85rGfhDVlcRPAOSHEhpTyMfvus7Cq8rP2308DOCOE2IZVvT/MOezholoxUZju6Ea1kqoNM/JtJaJp8JrHfhnAZQCfHXLfub6/l4YdR+HR2f0o4rNi1Fo31XoLWJzzyRCFENdjTxA1kyTqFbtqJXEAlWg4BnuCqLnfYV9SwIt6Y2KwEw3HYE8QVbFHfvC0Z9s9IhrEYE8Qo2FiIZNGOuW+32nYdSt2zoohGobBniCVeisWs0gKbMUQuWKwJ4jRMCO/TgwAFNTgaYMVO9EwDPYEqdbNTihGWb4z3ZEVO9EwDPYEqdZbsajYM+kUdC3FYCdywGBPkDhdqVnMapwVQ+SAwZ4gcWnFANaa8pwVQzQcgz1BjEYr8hcnKQUu3UvkiMGeIJV6fFoxBbZiiBwx2BNCSgmjEY/BU8C6+pStGKLhGOwJUTfbaLVlbCr2YpatGCInDPaE6O6eFI9gLzDYiRwx2BMiLvudKgU9zStPiRww2BNCDTRGfb9TRVXsUnLfU6J+DPaE6LRi4lKxZzWYbYlGqz3vUyEKHQZ7Qqi2RVxmxahPHpwZQzSIwZ4QcazYAS4ERjRMPF7l5CmOs2KAeOyi1Gy18bFf/X+4drvm+zECwCd/8F78yPfcNb0To8iKx6ucPFU72+LFpBUTo4p9p1zHU9s38d4jyzi6mvf1mC89fx1/+OIOg52GYrAnRLfHHo+nPE49drUX7cPft4Effs+dvh7zl/79H6Ecgzc1mo54vMoT7MpOBV+5ctPzuK+9tIt0SiCrxWNYJU4Ve7nWBAAUc/5fjotZrfM4on4M9oj7zO8/j4vPv+Xr2I31AoSI9kbWSmdD6xhcpFSuWW9OxRE+TRVzGt4q++/JU7Iw2CNub7+JzXtW8J/+zvs8j11aiM/Tndfjsz2eCvalUSr2nIYrO9H/t9N0xOeVnlCVegt3LS9gfTE771OZqTjNilE99lFaMcWshkot+v92mo54NFwTzGiYsbnoaBRZLYV0SsSiYlcBvZjL+H7MYi7DwVNyxGCPuGq9hXxM5qaPQghhLQQWg1kx5VoTQgD5jP836MWchobZRt2M/r+fgsdgjzijYcZmYa9RxWVN9nLdRFHXkEr5H9hWA61sx9AwySv1YqTdVrsiJfNpzGc1fOXKTfy9zz3Tue2j7zqMv/Ked8zxrEZXrplYHKG/DqBzfLlm4lAxWeMr5C2ZiRATRjNeC3uN6qPvOoz/+Y1reO71PQDAm3s1vHbLiFywV2rmSAOnQE/FHoNPLBQ8BnuEGZ1lApL5NP78R78DP//R7+j8/Wcfv4RvXa/M8YzGU6mbIw2cAt2B1jJbMTQEe+wRpi7OGeXCljgrZrVIBl251hz5Oey2Ynj1KQ1isEdY3Bb2mlQxm4lka6JcZyuGgsVgj7C4rbE+qcWchkrdRLsdre3yyjVzpKtOgYODp0T9GOwRZtitGFbsFhV2lYhdjVqpmSO3YlSFz4qdhmGwR5i6nJ49dksU53abrTb2m62RB0+zWhq6lsJt9thpCAZ7hHV67Ax2ANGsYjvrxIzxHC5yvRhywGCPMHU5fVKvPO3XnQIYnSq23FknZvRgL9pjCkT9GOwRZjSSPY+9n6p6ozSgOEmwL+aiOb2Tpo/BHmGVegt6OgU9JrsiTWox0q2Y0Xrs1mPYiqHhmAgRZjRM5BO6nMAwUZwCWKlbbaPxKvbMxIOnb+7t4/k3b0/0NSh8GOwRVq23UGAbpiOKs2I62+KNE+zZyXvs/+r3/wwf/7Wtib4GhQ+DPcKMhsk57D0KugYhEKkNKCbtsU8a7N+6Xsbrpf1YLH9MXQz2CKvUTV512iOVEijqWjRnxYzTY7cHT6Uc70rbdlvi6o0qAOClm9WxvgaFE4M9wqy12Fmx9yrmojWgWKk3kU4J5DKjvxSL2QxabYlasz3W937zdg1103qsCniKB89yTwixDOAMgG0AGwAuSikvOxz7KQCHAPx3AKsATkspzwZ2tnRAtW5itZCf92mESjGAvvMsVexNNoTwv3uS0rvC48IYLbmrO90w395hsMeJn8/x5wGclVJuA4AQ4oIQ4rSUsuRw/Bn7z0UADwdyljSU0Wjx4qQ+xYjN7S6PsU6M0gn2uok7xnj89g1r7fq8nmbFHjOuv1F2tb6hQt22DeAUgCeHPKQkpVwJ7vTITZU99gGLuQz29iPUYx9jkw2lM29/zDey7Z0qCnoa7z26jG0Ge6x4NfY2AZT6bisBeNDtQUKIk0KIjfFPi/yoNhjs/az1UyIU7LUmFsd8DtVFTeN+Qrl6o4rj6wWcWC9ie6cy9iAshY9XsC8D2O277Sas/vlQQoiHYFX1J4UQj7gcd0YIsSWE2NrZ2fF5uqSoQTNOdzwoarsoVcbYZEPpbrYx3hvZ1RtVHF8r4vhaAeWaiZvVxlhfh8LHz1C8Y4j3k1I+JqV8UkpZklI+CeAhIcQpl2M3pZSb6+vrfr8F2Qwu2TtUEHO7Z0kNno5DPe72GG9kdbOF124ZOL5WwPG1AgDOjIkTr2Avwaraex3CYBUPwGrB9N10GR5tGxqPWtmRC4AdVMxpMBottCKyi1IQg6fj9NhfuWmgLYET6wVsrBUBHJwlQ9HmFexbGKzYlwFc6D/QDvUvDTn2ypjnRi7UJhucx35Q1JYVmGTwdJJ9T9Vg6fG1Au5aWYCeTuGKPUuGos812O0pjVt9A6GbsKYyQgixoe6z57af6/sSGwCeCOxsqcNgxT7UklqTfcy+8yzVzRYaZnvsVoyWTmEhkx7rSturPcGeTgnccyjPij1G/PxGnQZwRgixDat6f7hnDvtZWFW5ughpy75IqQTgBKwLlEqgwFXqrNiHidIuSupTxSTjJMWchvOXXsMffevGSI97q1zH+mK282nh+FqBUx5jxPM3yg7mzzrcd67v75dh9dVpytTgKVd3PChKm21Msi2e8jMfOYGnt2+O/Lijq3l8+N61zt/vWMriay8NHTqjCGIqRFS1YW+Lx4r9gOKEF+3MUvdT1/gvw5/+8HH89IePT3wuhazWGZCn6OMiYBFlqI2sWbEfsNRzmX3YGfabcximrBZ1DY1WGw1zvAXFKFwY7BEVRLUXR92rMcM/eKqewzDsgqV+j1SLj6KNwR5RqtrjlacHRakVUw2gxx6USaZOUvgw2COq2jChaylk0nwKexX0NISIRkBVQ/SpS50D++zxwFSIKKPOJXuHEUJEZr2Yih2ixRCMk6hB+Ci8IZI3BntEccleZ4sRCfZqiK5F6Fbs4f+5kTcmQ8RIKVE327hdMzmH3cFiLjP2ioezVK2byGopaCFop6nfJQ6exgOTIWJ+7oln8T+eeR0AsHkP9zQZppjT8NINA5//+usD9913eBHfeefSHM5qUJg2I+8OnrLHHgfh+K0i3164XsY77yjir528Gw+cODTv0wmlu5YX8DvPvoG//xtfH3rfH/+jH5j9SQ0Rps3I1XmwFRMPDPaIMRotfPddb8Mnvv/EvE8ltP7t6T+Hf3Dq3oHbf/2rr+A///FV1Jot5DLzD9RKPTzttAKnO8ZKOH6ryLdq3eRsGA+6lsLGenHg9nffbbVgXi/t48SQ+2etWh9/LfagZbUUtJRgxR4T8x+1oZFYH9/DEQZRc2QlDwB4ZdeY85lYwjSzSQhhrxfDYI8DBnuESCmtDaxZsY/l6KoV7K+FJNgrIarYAeviLg6exgODPUL2my1ICeRDFAZRsr6YRVZLhahiD8/gKWD12TndMR4Y7BGiLvdmxT4eIQSOrOZDFOxmqFbnLGSjtRE4OWOwR4iqpsIUBlFzZGUBr+7uz/s0Om21MLViiuyxxwaDPUK4VO/kjq7m8equASnlXM+j1myjLcP1XBayaS4CFhMM9ggxuGvSxI6s5lGum9jbn++SA91t8cLzXLIVEx8M9gipctekiR1ZDceUxzAt2asUdA1VDp7GAoM9QlixT07NZZ93nz2MbbVCVoPBVkwsMNgjpFPlsWIf25HVBQAhqthD9FwWs2nuexoTDPYICePH96hZzGWwks/g1VtzDvZGeNZiV7gme3ww2COkyn1OA6FmxsxTZ/ekEL1JcyGw+AjPbxV5Mhom0imBrMb340ncvZrH/37uGj78yJc7t73rHUt49Cc2Z3YORgg/fak3GQ6gRl94fqvIU9Xe51QIMe9TibSf/OAx5LQ0JKy57C9cK+OLz12H2WrPbDejsA6eAmzFxEF4fqvIU5hWA4yyzWOr2Dy22vn7f3nqJfzTzz+H0n4Ta8XsTM4hjMtDqHPhQmDRx8/0EWI0WuyvT8FyXgcAlIzGzL5ntRGe/U4VVTQYrNgjLzy/VeSp2mDFPg2rdrDvVmd3NWrYluwFevc9ZbBHHYM9Qow6K/ZpWClkAAC71RlW7CFsq7HHHh8M9ggJY5UXB6sFq2K/NctWTCiD3d7QusEee9Qx2CPEaIRr/e64WMnPPtitN+lwffrKamlk0oKtmBhgsEdItRGuHXfiIpdJYyGTxq0ZtmLCunct9z2NBwZ7hBgh23EnTlYL+swHT0MZ7LrGNdljgMEeEe22hNEMZ5UXB8v5zOx77CEcCLc222DFHnUM9ohQG1mHMQziYLWgzzjYw/kmXchyTfY4YLBHhHqx5UMYBnGwktdn1mMP436nSpG7KMUCgz0ijBBegh4nVo99NsFuNOxPXyEMdqvHzmCPuvD9ZtFQ3fW7+ZRNw3I+g9s109dCYEbDxO/+yZtotsbbkCLM6+oXshp2ynU8/vTLjse8Y3kBf/H+O2Z4VjSq8P1m0VDdRaP4lE2DukjJz0Jgv/3MG/jHv/WNib/nUXv/1TA5diiPW0YT/+S3vul4jBDA5V94ECv2z4zChykREd0eO1sx09C5SKna8Az2P7t2G8Wshi/93Ecw7gLKmXQqlMH4d3/gnfjxDxyBvaLxgK9e3cUnP/cMXrhexp/fODTbkyPfGOwRYbBinypVsfvps79wrYz7DhdxeCk37dOaOSEE7lh0/ne9/9gKAOBFBnuocfA0IjoVOwdPp2I5by0E5jXlUUqJF6+Xcf/bF2dxWqHz9qUcFnMaXrxenvepkAsGe0SoAbcwTpGLg+5CYO5Xn96oNHDLaOLeO5IZ7EII3Hd4ES9eq8z7VMgFgz0iDLWRNXvsU7GS99eKUZVqUit2ALjv8CJeuF6GlA6NeJo7BntEVOsmtJSAHqIdd+Ikl0kjr3svBPbCNSvY7zuc3GC//3ARe/tNvFWuz/tUyAFTIiLUtnjcyHp6VvI6dj167N96q4zVgo61YvhmtMzKffanFfbZw8uzYSuEWAZwBsA2gA0AF6WUlyc9lkZT5SYbU7dSyKDk0WNXM2KS/AZ7v/1p5YVrZXzfvetzPhsaxk9SnAdwVkq5DQBCiAtCiNNSytKEx9IIqg2T68RM2UrefVkBa0ZMBT928q4ZnlX4HCpmcaigs2IPMddWjF2Bb6igtm0DODXJsTS6ar3FdWKmzGuFxzf2aqjUzUT31xVrAJUzY8LKqwTcBFDqu60E4EEAT05wbCA+//XX8c3X96bxpUPnxetlHDtUmPdpxNpKXse1vRp+6ff+dOj9b+7VACR7Roxy/9sX8Rtfe8XxZ0X+feyBYzgS8PISXsG+DGC377absPrnkxwLIcQZWP14HD161OM0hvvKt2/iC3/yxliPjaL33bMy71OItZP3rODJS6/h8adfcTzm6Goe33nn0gzPKpw+ct86fvOy+8+K/Pmh774z8GAXbnNRhRAPAfi0lPJ9Pbd9CsD7pZSnxz223+bmptza2hrzn0BElExCiEtSys3+272mO5ZgVeK9DmGwMh/1WCIimhKvYN8CsNp32zKACxMeS0REU+Ia7PY0xS0hRG+ffBPARQAQQmyo+7yOJSKi2fAzMfo0gDNCiG1YFfnDPfPSz8Kqys/6OJaIiGbAM9jtYP6sw33n/B5LRESzwbViiIhihsFORBQzDHYiophxvUBpZichxA6Al8d8+BqAGwGeDgWDz0v48DkJn0mfk3uklANLbIYi2CchhNgaduUVzRefl/DhcxI+03pO2IohIooZBjsRUczEIdgfm/cJ0FB8XsKHz0n4TOU5iXyPnYiIDopDxU5ERD0Y7EREMRPa3ZGFECdhrQ65DWuhsZKU0nOlyHEfR/6M8/O1V/w8D+BRWMs7b8DaMvFRKeXlqZ5wQgghzgP4jN+fJ18nszHK8xLk6ySUwW7/Az/du/OSEOK8EGLX7R847uPInwl/vssAHrH/exHAOT4nk7Gfj3OwNrk5BSsQ/D6Or5MpGfd5sS0jiNeJlDJ0f+wfxKm+204CuDCNx/HP1J+XDQAn533+cf4D4Er/cxP088g/U39eAnudhLXH/tdhfUTstQ3r3W8ajyN/+PONBz6PMRe6YLc/xixLKQ/84kl7ww67NxjY48ifIH6+QohlIcTJvl22aIb4Ogm/IF4noQt2DG6I3a9/X9VJH0f+LHvc7/Xz/XH0DNbZPV2vr0nBW/a4n6+T+QrkdRLGYKeYsavDz0gpL0opS9IaDLoA4FfmfGpEoRHk64TBTjMhB/e+vQjgIVbtRF1BvU7CGOwlwOozOdy/G/DjyJ8SMN7PVwhxxuV49ttnqwTwdRJGQb5OQhfs9seREvp6feoXUTrM6Rz3ceTPuD9fewDoUZeBoP7ZGTRFfJ2EU9Cvk9AFu+0irHm1vTbt26fxOPJn5J+vHSRn+2dhwJpyd3nIR0+aPr5OQibo10lYg/0cgE/33XbWvh1AZ0rQJSHEQ6M8jiYy7vOy21uJ2NXhWQAPT/FcCXydhNW0XyehXbZXCHEKVlWxDau/dFn2rGVh/6Ovwrrk9jG/j6PJTPC8PIRun/AEgEeGVCc0Avtn/WlYUxjPALgMq+q+oJ4Tvk5mb8LnJZDXSWiDnYiIxhPWVgwREY2JwU5EFDMMdiKimGGwExHFDIOdiChmGOxERDHDYCciihkGOxFRzDDYiYhihsFORBQz/x+hqdm6j5sYmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xseq = np.array(range(0,151))/100\n",
    "plt.plot(xseq, es_normal[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ES function is usually quite similar to the Betti curve although it can be completely different for others. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stability of the ES-function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The main theoretical advantage of the ES function respect to the Betti curve is that it is more robust to modifications in the input. Actually, the following result holds:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ |ES_1 - ES_2| \\le d_\\infty(B_1, B_2)\\left(\\log(n_{\\infty})-\\dfrac{3}{2}\\log\\left(\\dfrac{3}{2}\\dfrac{n_\\infty}{L_\\max}d_\\infty( B_1,B_2 )\\right)\\right) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $ES_1$ and $ES_2$ are the ES function of the barcodes $B_1$ and $B_2$, $L_\\max$ is the maximum total length of the bars in the barcodes, $d_\\infty$ is the bottleneck distance and $n_\\infty$ is the sum of the numbers of unpaired bars and the number of pairs given by $d_\\infty$. This formula can be directly deduce from the one appearing [here](https://arxiv.org/abs/1803.08304)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you might notice, the source of unstability might be $L_\\max$ and $n_\\infty$. Nevertheless, in opposition to what happened with PE, these terms are inside a logarithm, making it much more robust to noise. This is an advantage with respect to, for example, the Betti curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When should I use the ES function?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually, the ES function can be used in the same fashion than the Betti curve. If you expect noise to introduce artificial small bars in the barcodes, then you should use the ES function instead. When the barcodes have some bars which are considerably bigger in comparison with the others, the ES function and the Betti curve tends to produce different values. In that cases, both curves can complement each other. Please, refer to this [paper](https://arxiv.org/abs/1904.07768) to a more general study of persistence curves, including the Betti curve and the ES function (denoted there as life entropy)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall you can use the ES function as a parameter in your scikit-learn pipeline, together with other methods from the [representations](https://gudhi.inria.fr/python/3.1.0.rc1/representations.html) module."
   ]
  }
 ],
 "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}