{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison to SOTA\n",
    "\n",
    "This notebook contains a comparison of the Vietoris-Rips (VR) persistent homology (PH) computation (on the GPU) implemented in `torchph` to the implementations available via `dionysus` (v2) and `ripser` (which is SOTA in terms of runtime).\n",
    "\n",
    "This runtime study is a small-scale version of the runtime study in \n",
    "\n",
    "**Connectivity-Optimized Representation Learning via Persistent Homology**    \n",
    "C. Hofer, R. Kwitt, M. Dixit and M. Niethammer    \n",
    "ICML '19    \n",
    "[Online](http://proceedings.mlr.press/v97/hofer19a.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: ripser in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (0.4.1)\r\n",
      "Requirement already satisfied: numpy in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from ripser) (1.18.1)\r\n",
      "Requirement already satisfied: scipy in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from ripser) (1.4.1)\r\n",
      "Requirement already satisfied: persim in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from ripser) (0.1.2)\r\n",
      "Requirement already satisfied: Cython in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from ripser) (0.29.15)\r\n",
      "Requirement already satisfied: scikit-learn in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from ripser) (0.22.1)\r\n",
      "Requirement already satisfied: matplotlib in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from persim->ripser) (3.1.3)\r\n",
      "Requirement already satisfied: hopcroftkarp in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from persim->ripser) (1.2.5)\r\n",
      "Requirement already satisfied: joblib>=0.11 in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from scikit-learn->ripser) (0.14.1)\r\n",
      "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from matplotlib->persim->ripser) (2.4.6)\r\n",
      "Requirement already satisfied: python-dateutil>=2.1 in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from matplotlib->persim->ripser) (2.8.1)\r\n",
      "Requirement already satisfied: kiwisolver>=1.0.1 in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from matplotlib->persim->ripser) (1.1.0)\r\n",
      "Requirement already satisfied: cycler>=0.10 in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from matplotlib->persim->ripser) (0.10.0)\r\n",
      "Requirement already satisfied: six>=1.5 in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from python-dateutil>=2.1->matplotlib->persim->ripser) (1.14.0)\r\n",
      "Requirement already satisfied: setuptools in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib->persim->ripser) (45.2.0.post20200210)\r\n"
     ]
    }
   ],
   "source": [
    "# LOCALLY INSTALL RIPSER\n",
    "!{sys.executable} -m pip install ripser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: dionysus in /scratch4/chofer/opt/anaconda/envs/pyt_1.4/lib/python3.7/site-packages (2.0.6)\r\n"
     ]
    }
   ],
   "source": [
    "# LOCALLY INSTALL DIONYSUS\n",
    "!{sys.executable} -m pip install dionysus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "from collections import defaultdict\n",
    "\n",
    "# IMPORT RIPSER and DIONYSUS\n",
    "from ripser import ripser\n",
    "import dionysus as d\n",
    "\n",
    "import torch\n",
    "import time\n",
    "\n",
    "from scipy.spatial.distance import pdist\n",
    "\n",
    "from sklearn import datasets\n",
    "\n",
    "from torchph.pershom import vr_persistence_l1, vr_persistence\n",
    "\n",
    "device = 'cuda'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluating N=50\n",
      "Evaluating N=100\n",
      "Evaluating N=150\n",
      "Evaluating N=200\n",
      "Evaluating N=250\n",
      "Evaluating N=300\n",
      "Evaluating N=350\n",
      "Evaluating N=400\n",
      "Evaluating N=450\n"
     ]
    }
   ],
   "source": [
    "def gen_random_10D_data(n):\n",
    "    return np.random.randn(n,10)  \n",
    "\n",
    "def gen_circlces(n):\n",
    "    return datasets.make_circles(\n",
    "        n_samples=n, \n",
    "        noise=0.05, \n",
    "        factor=0.5)[0]\n",
    "\n",
    "N = np.arange(50,500,50)\n",
    "\n",
    "times = defaultdict(list)\n",
    "\n",
    "f = gen_random_10D_data\n",
    "\n",
    "for n in N:\n",
    "    trials = defaultdict(list)\n",
    "    print('Evaluating N={}'.format(n))\n",
    "    \n",
    "    for i in np.arange(50):\n",
    "        \n",
    "        #Generate data\n",
    "        x = f(n)\n",
    "        D_l1 = pdist(x, metric='cityblock')\n",
    "        thr_l1 = np.max(D_l1.ravel())\n",
    "\n",
    "        X = torch.Tensor(x).to(device)\n",
    "        \n",
    "        # Timing torchph\n",
    "        t0 = time.time()\n",
    "        l, _ = vr_persistence_l1(X.contiguous(),0, 0);\n",
    "        t1 = time.time()\n",
    "        trials['torchph'].append(t1-t0)\n",
    "        \n",
    "        if 'ripser' in sys.modules:        \n",
    "            t0 = time.time()\n",
    "            ripser(x, maxdim=0, thresh=thr_l1, metric='manhattan')\n",
    "            t1 = time.time()\n",
    "            trials['ripser'].append(t1-t0)\n",
    "        \n",
    "        if 'dionysus' in sys.modules: \n",
    "            t0 = time.time()\n",
    "            dd = pdist(x, metric='cityblock')\n",
    "            filt = d.fill_rips(dd, 1, thr_l1)\n",
    "            m = d.cohomology_persistence(filt)\n",
    "            t1 = time.time()\n",
    "            trials['dionysus'].append(t1-t0)\n",
    "        \n",
    "    times[n] = trials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "methods = []\n",
    "if 'ripser' in sys.modules:\n",
    "    methods.append('ripser')\n",
    "if 'dionysus' in sys.modules:\n",
    "    methods.append('dionysus')\n",
    "methods.append('torchph')\n",
    "\n",
    "T = np.zeros((len(N),3))\n",
    "S = np.zeros((len(N),3))\n",
    "for j,(k,v) in enumerate(times.items()):\n",
    "    for i, method in enumerate(methods):\n",
    "        T[j,i] = np.array(v[method]).mean()\n",
    "        S[j,i] = np.array(v[method]).std()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEYCAYAAABPzsEfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUVfrA8e+bmfRGCCSQQgIqqCC9qUhREURFRSy4YkNZXFkLPztrd9ddde3Y+8qKioqssqKwgliQJopUqRIghBpIn3J+f9wZMgkTMkMyqe/nee6TW86deSeB886559xzxRiDUkopVZ2w+g5AKaVU46AJQymlVEA0YSillAqIJgyllFIB0YShlFIqIJowlFJKBUQThlJKqYBowlBKKRUQTRhKqVonIh1E5HURme6zL1ZE3haRV0XkD/UZnzo6mjBUnRMRl4gsF5GVIvKziEwSkSr/LYpICxH5Uw3fM1tEfq207wERuS2Ac7+vLgYRSRWRf4vIRhFZKiI/iMiFNYk5UCLyfUN7LWPMRmPMuEq7RwHTjTHXAyNr431U3dKEoepDsTGmuzGmMzAUGAHcf4TyLYAaJYyaMMaccqQYRESAGcA3xpgOxphewGVARh3GVy+vJSJtReQNEXlLRF4TkRdF5LgqimcAWz3rrhoFquqFJgxVr4wxecB4YKJYZni+oa8UkfGeYn8HjvG0Sh4HqKLcUfG0PlZ7LpWsFJEvRSTa53iBvxh8nA6UGWNe8vlcW4wxz/m8xmHxVm71iMhtIvKAZz1WRD73tMB+FZFL/e3zia+696ny81X6XRQEUx4YArwBvGSMuQ74Bji7irI5lCdRrXsaIf2jqXpnjNmI9W8xBbjW8w29N3CTiCQDdwEbPK2S2z2n+SuHiMwSkbSjCOM4YIqn1bMfuKjScX8xeHUGllXz+n7jPYLhwHZjTDdjTBfgiyr2Bfo+1X2+ygItHwMcDyz2bIcBBSKSLCIvAT1E5G7PsY+Bi0TkReA/1by/aoDs9R2AUh7i+XmTz7X/TKyKK9dPeX/l9hhjRlTx+lVNy+zdv8kYs9yzvhTIDjTwykRkCjAAq9XR5wjx+vtcXiuAJ0TkH8BnxpgFInLYPj/nVfU+wX6+QMuXAbHGGJeIxAH9gBXGmD3ABN+CxphC4Jpq3lc1YNrCUPVORDpgXdM+ETgTONkY0w34CYjyU35wIOUq2QMkVdrXEtjtWS/12e8iuC9TK4Ge3g1jzI3AGUDrauJ1UvH/YJTPa6wDemEljkdF5D5/+3yDqOb3Euznq7a8iBwPbAZGeTrLPwEWG2Nerea1VSOlCUPVKxFpDbwEPA8kAPuMMUWeyqi/p9hBIN7ntMQqylXJGFMA7BCRMzzv2xLrEs+3AYZaOQZf/wOiROQGn30xAcS7E0jxXL6JBM71nuC5rFZkjHkXeALo6W9fpTiC/r3UUGespHuzZ3kG6/ekmii9JKXqQ7SILAfCsb5l/wt40rM9QUR+AdYCCwGMMXtE5DtPB/F/gb/4KwdWHwZwnTFmu5/3vRKYIiL/9Gw/aIzZICLZ1QVcOQbffgxjjBGRC4CnROQOYBdQCNzpKfJFFZ/LISIPAT8Cm4A1Pm95EvC4iLgBB3BDFft8+X2fEBoIZBhjLgIQkaeB40XkW2PM7iOfqhoj0SfuKaWUCoReklJKKRUQTRhKKaUCoglDKaVUQDRhKKWUCogmDKWUUgFpssNqW7VqZbKzs4/6/MLCQmJjY2svoFqicQVH4wqOxhWcphjX0qVLdxtjWvs9aIxpkkuvXr1MTXz99dc1Oj9UNK7gaFzB0biC0xTjApaYKupVvSSllFIqIHWaMERkuIisFZH1InKXn+MTRGSFZwrpb0XkRJ9jd3vOWysiw+oybqWUUnWYMETEBkzBmiv/RGCMb0Lw+Lcx5iRjTHfgMazpIvCUuwxr7prhwAue11NKKVVH6rLTuy+w3ljPPkBEpgHnA6u8BYwxB3zKx1I+9fT5wDRjTCmwSUTWe17vh2ACcDgc5OTkUFJSUm3ZxMREVq9eHczL14n6jisqKoqMjAzCw8PrLQalVP2oy4SRTvnjGcF6+la/yoVE5EZgEhCB9SQz77m+E6nlePZVPnc81tPbSE1NZd68eRWOx8XFkZqaSnp6OiJS+fQKXC4XNlvDa8TUZ1zGGPLz8/n5558pKCiocKygoOCw33dDoHEFR+MKTrOLq6re8NpegIuB13y2xwLPHaH85cDbnvUpwBU+x14HLjrS+/kbJbVq1SrjdrsDGilw4MCBgMrVtfqOy+12m1WrVh22vymOFgkljSs4GldwmsIoqRysJ4B5ZQD+pqD2mgZccJTnVqm6loU6Mv39KdV81WXCWAwcJyLtRSQCqxN7pm8BETnOZ/Mc4DfP+kzgMhGJFJH2WI+dXFQHMSullPKos4RhjHECE4HZwGrgA2PMShF5SERGeopNFJGVnofrTAKu8py7EvgAq4P8C+BGY4yrrmKvTTabje7du9OlSxfOO+889u/fD8D27dsZPXp0PUenlGrsjDG4QlQ91ul9GMaYWcaYjsaYY4wxf/Xsu88YM9OzfrMxprMxprsxZognUXjP/avnvE7GmP/WZdy1KTo6muXLl/Prr7/SsmVLpkyZAkBaWhrTp08P2fs6nc6QvbZSquHYV7IPp9vp7e+tVXqndz06+eST2bZtGwCbN2+mS5cuALz11lucf/75DB8+nE6dOvHggw8C1vwwo0ePplu3bnTp0oX3338fgKVLlzJo0CB69erFsGHD2LFjBwCDBw/mnnvuYdCgQTzzzDP18AmVUnWpxFnC/M3zySnKCcnrN9nJB6t1yy2wfHmVh6NdLgh2+Gr37vD00wEVdblczJ07l3Hjxvk9vmjRIn799VdiYmLo06cP55xzDlu2bKFt27bMnj0bgPz8fBwOB3/+85/59NNPad26Ne+//z6TJ0/mjTfeAGD//v3Mnz8/uM+hlGp03MbN+j3rmTR7Eu4yN1efczW2Wr6/WVsYday4uJju3buTnJzM3r17GTp0qN9yQ4cOJTk5mejoaEaNGsW3337LSSedxLx587jzzjtZsGABiYmJrF27ll9//ZWhQ4fSvXt3HnnkEXJyyr9dXHrppXX10ZRS9Whv0V4e/OZBtuRvYeKxEwmT2q/em28Lo5qWQPHBg8THx9f623r7MPLz8zn33HOZMmUKN91002HlKg9fFRE6duzI/PnzWbBgAXfffTdnnXUWF154IZ07d+aHH/zf9N4Qp15WStWuYkcx7/36HtNXTWdin4l0jewakvfRFkY9SUxM5Nlnn+WJJ57A4XAcdvyrr75i7969FBcXM2PGDE499VS2b99OTEwMV1xxBbfddhvLli2jU6dO7Nq161DCcDgcrFy58rDXU0o1TS63i+U7lnPfvPs4KeUk/u+U/wvZezXfFkYD0KNHD7p168a0adM47bTTKhwbMGAAY8eOZf369Vx++eX07t2b2bNn83//93/Y7XbCw8N58cUXiYiIYPr06dx0003k5+fjdDq55ZZb6Ny5cz19KqVUXcorzOP2ObdT4izh+RHP4zZuBAnJTbaaMOpY5TmY/vOf/xxa//XXXw+tp6Sk8Pzzz1coO2zYME455ZTDLpV1796db7755rD3aohz3Cilak9hWSHP/vgs3239jkfPeJT2LdpT5Cgi3BaayUH1kpRSSjVCTreTeZvm8dTCpxjaYShju46l0FFI27i2CKGZwkdbGA3Q1VdfzdVXX13fYSilGrCt+VuZ9NUkEiITeOKsJyh2FJMYmUhCVELI3lMThlJKNTIHSg7wwPwHWLdnHf+68F8kRCbgdrtJiU0J6fvqJSmllGpEHC4H01dN552f3+Ga7tcwOHswJc4S0hLSsIWF9lk5mjCUUqqRMMawatcq7pp7Fx2TOzL5tMkUlhWSGptKlD0q5O+vl6SUUqqR2F+yn9u+vI39Jft5d9S7AMSEx9AiqkWdvL+2MOqYd3rzzp07061bN5588kncbjcAS5Ys8XvXt1JKlTpLmbJ4CnM2zeGuAXdxfKvjcRs3bePb1tmDzbSFUce8U4MA5OXlcfnll5Ofn8+DDz5I79696d27dz1HqJRqaNzGzQ9bf+BvC/7GqZmncn3P6yl0FJKZkIk9rO6qcW1h1KOUlBReeeUVnn/+eYwxzJs3j3PPPReAvXv3csEFF9C1a1f69+/PL7/8AsDf/vY3rr32WgYPHkyHDh149tlnAbj33nsrTGE+efJknn32WXbs2MHAgQMPPbRpwYIFAMTFxR0qO3369EPDeD/88EO6dOlCt27dGDhwYF38GpRS1dhVuIubvriJSFskTw9/mhJnCcnRycRG1O1ccc22hXHLF7ewPLfq6c1dLhe2IKc3796mO08PD2x6c68OHTrgdrvJy8ursP/++++nR48ezJgxg//9739ceeWVh1oma9as4euvv+bgwYN06tSJG264gXHjxjFq1Chuvvlm3G4306ZNY9GiRbz11lsMGzaMyZMn43K5KCoqOmI8Dz30ELNnzyY9Pf3Q0wCVUvWnxFnCg/MfZEXeCl469yVaRrckjDCSY5LrPJZmmzAaEn9Pxvr222/56KOPADj99NPZs2cP+fn5AJxzzjlERkYSGRlJSkoKO3fuJDs7m+TkZH766Sd27txJjx49SE5Opk+fPlx77bU4HA4uuOACunfvfsRYTj31VK6++mouueQSRo0aVfsfVikVMLdxM3PtTF5e+jKXdL6EEceOoMhRRPuk9iGZvrw6zTZhVNcSOBii6c0r27hxIzabjZSUFFavXn1ov78k4u3YioyMPLTPZrMdevzqddddx1tvvUVubi7XXnstAAMHDuSbb77h888/Z+zYsdx+++1ceeWVFTrJSkpKDq2/9NJL/Pjjj3z++ed0796d5cuXk5xc999klFKwad8mJs2eRGZCJg8PeZhCRyFpcWlE2CLqJR7tw6hHu3btYsKECUycOPGwUQ4DBw5k6tSpgDWJYKtWrUhIOPIt/xdeeCFffPEFixcvZtiwYQBs2bKFlJQUrr/+esaNG8eyZcsASE1NZfXq1bjdbj755JNDr7Fhwwb69evHQw89RKtWrdi6dWttfmSlVICKHEVM+nISuQW5PHv2s4QRFvKpP6rTbFsY9cX7xD2Hw4Hdbmfs2LFMmjTpsHIPPPAA11xzDV27diUmJoa333672teOiIhgyJAhtGjR4lD/y7x583j88ccJDw8nLi6Od955B4C///3vnHvuuWRmZtKlS5dDs+jefvvt/PbbbxhjOOOMM+jWrVstfnqlVCBcbhevLX2NmWtnMqn/JLqmdq2TqT+qowmjjrlcriqPDR48mMGDBwPQsmVLPv3008PK3HPPPRUulflOie52u1m4cCEffvjhoX1XXXUVV1111WGvM3r0aEaPHn3Y/o8//jigz6GUCp1lO5Zx77x76dm2Jzf1u4liZzFZiVkhn/qjOnpJqolYtWoVxx57LGeccQbHHXdcfYejlDpK+4v3c+OsG3G5XTx39nOUOEtIiU0hOjy6vkPTFkZTceKJJ7Jx48b6DkMpVQNOt5O/Lvgri7cv5smzniQ1NpUIWwRJUUn1HRqgLQyllGoQjDHM2TCHZ358hhHHjWDUCaNwGzdt4trU2dQf1anThCEiw0VkrYisF5G7/ByfJCKrROQXEZkrIlk+x1wistyzzKzLuJVSKtR2FuzkT7P+RHJMMv848x8UO4tJi08L2eNWj0adXZISERswBRgK5ACLRWSmMWaVT7GfgN7GmCIRuQF4DLjUc6zYGHPku86UUqoRKnOVMenLSWzav4lpo6cRaYukRVSLOp/6ozp12cLoC6w3xmw0xpQB04DzfQsYY742xnjnrlgIZNRhfEopVeeMMUz9ZSrv/foef+z1R/qm9SU8LLxepv6ojvi7ozgkbyQyGhhujLnOsz0W6GeMmVhF+eeBXGPMI55tJ7AccAJ/N8bM8HPOeGA8QGpqaq9p06ZVOJ6YmMixxx4bULxHM5dUdfbs2cPIkSMB2LlzJzabjVatWgHw9ddfExFR/d2b/uJ66KGHSE5O5sYbbww4luOPP56FCxfSokXw8+ivX7/+0DQlXgUFBRUmNGwoNK7gaFzBqY24dpfu5rql15Eckcyz3Z/Fjp0IewTC0fdb1CSuIUOGLDXG+J02uy5HSfn79H6zlYhcAfQGBvnsbmeM2S4iHYD/icgKY8yGCi9mzCvAKwC9e/c23nsavFavXh3wdB+hmBokPj7+0KyzDzzwAHFxcdx2220BnWuMwRhDYWHhYXFFRkYSFRUVVLwiQnx8/FF9xqioKHr06FFh37x586j8+24INK7gaFzBqWlcJY4Shk8dTom7hNcueo228W1Ji0ur8d3cofp91eUlqRwg02c7A9heuZCInAlMBkYaY0q9+40x2z0/NwLzgB6Vz23MHnvsMbp06UKXLl147rnnAOubfJcuXZgwYQI9e/Zkx44dfPHFF/Ts2ZNu3bpx1llnHTp/xYoVDBo0iA4dOjBlypRD53fu3JmxY8dy0kkncckll1BcXHzonKeffpoePXrQtWtX1q1bV7cfWKlmzm3cPP7948zfMp+/DPwLmQmZ9T71R3XqsoWxGDhORNoD24DLgMt9C4hID+BlrEtXeT77k4AiY0ypiLQCTsXqED9qt9wCy6ue3RyXK5pgr0h17w5PBze7OQCLFi1i6tSpLFq0CJfLRd++fRk0aBAxMTGsWrWKN998k5deeonc3FwmTZrEt99+S1ZWFnv37j30GuvWrWPu3Lns37+fE044gQkTJgDWDX2vv/46/fv358orr+Tll1/mlltuAaz5pH766SeeffZZnnzySV566aXgg1dKHZWFWxfy1wV/ZUj2EK446Qrcpv6n/qhOnbUwjDFOYCIwG1gNfGCMWSkiD4nISE+xx4E44MNKw2dPAJaIyM/A11h9GKtoIhYsWMBFF11ETEwM8fHxXHDBBXz77bcAHHPMMfTp0weAH374gdNOO42sLGu0ccuWLQ+9xrnnnktERAQpKSm0bNmSXbt2AdC+fXv69+8PwBVXXHHodYFD05f36tWLzZs3h/xzKqUs+4v3c91/riM2IpZ/nvVPSl2ltI1vW+9Tf1SnTu/0NsbMAmZV2nefz/qZVZz3PXBSbcZSXUvg4MHiOpneHPxPZe4VGxtboVxVN/BUNeV55fK+295zfMsrpULL5XZxx5w7WL17NW+e/yYx4TG0jm3dIKb+qI7e6d0ADBw4kE8++YTi4mIKCgr49NNPOe200w4rd+qpp/LNN9+wZcsWgAqXpKqyadMmFi9eDMB7773HgAEDajd4pVRQZqyZwWvLXuOKrldwWrvTiAmPaTBTf1RHE0YD0LdvX8aMGUOfPn3o378/N9xwAyeddHiDKjU1lSeffJLzzz+fbt268Yc//KHa1+7cuTOvvvoqXbt2pbCwkPHjx4fiIyilApCTn8PE/06kfVJ7Jp82ucFN/VEdnXywnjzwwAMVtu+44w7uuOOOCvuOPfbYQ8/x9ho+fDgXX3xxhX2PPPJIhe01a9YA1igpm83GK6+8ctj75+TkHFrv378/c+bMCfozKKUC53A5+ONnf2R30W5mXmZ1zza0qT+qoy0MpZSqAy8ueZFZ62dx2ym3cUzSMSRHJze4qT+qoy2MJsxfC0UpVfdW7FzBPXPvoX96f8Z1H4c9zN4gp/6oTrNrYdTVVChNlf7+lApOsaOYq2dcjU1sPDX8Kdy4aRvfljBpfNVv44u4BqKiotizZ49WekfJGMOePXuIioqq71CUahSMMdz39X0sy13G3874Gy2iWtAmtg2R9sjqT26AmtUlqYyMDHJycg7d1HYkJSUlDbJirO+4oqKiyMjQSYSVCsT/Nv2PpxY+xYXHX8iwY4YRGxHboKf+qE6zShjh4eG0b98+oLLz5s07bIK9hqChxqWUqmhv8V7GzRxHm7g2PDD4AUSkwU/9UZ1mlTCUUqouGGP40+d/YuuBrXw4+kMibZGkxac1+Kk/qhNwH4aI3B/KQJRSqql455d3eH/l+9zY50Y6p3RuNFN/VCeYFsb9IhIDtASWAdOMMftCE5ZSSjVOG/Zu4Ob/3ky31G78qfefGtXUH9UJZpSUAUqwZpvNBL4XkW4hiUoppRohp9vJVTOuosxVxtPDn8YWZmtUU39UJ5gWxhpjjPey1HQReQt4CTi91qNSSqlG6NEFj/Ld1u/4+xl/p01cm0Y39Ud1gmlh7BaRXt4NY8w6oHXth6SUUo3Pom2LePibhxl2zDAu6HRBo5z6ozrBtDBuAqaJyFJgBdAV2BSSqJRSqhEpLCtk7CdjaRHVgkeGPEKkPbJRTv1RnYAThjHmZxHpDpwJdMF68t17oQpMKaUai1tn38q6Pet458J3SIhKaLRTf1QnqPswjDGlwOeeRSmlmr2Za2by6rJXubb7tfRJ69Oop/6oTkAJQ0T6AsYYs1hETgSGY3WCz6rmVKWUarLyHflc/5/r6ZjckVv730piZCLxkXXzaOf6UG3C8NywdzZgF5GvgH7APOAuEelhjPlraENUSqmGx+1289jax9hfsp93LniH6PBoUmJTmswQWn8CaWGMBroDkUAukGGMOSAijwM/ApowlFLNztM/Ps33e77nvoH3kd0iu0lM/VGdQHplnMYYlzGmCNhgjDkAYIwpBtwhjU4ppRqgT9d8yl1z7qJni55c1uUyUuJSmsTUH9UJJGGUeaYEATh0H4aIJKIJQynVzHy/9Xuu+OQKOiR14O7j7yYuIq7JTP1RnUASxkBP6wJjjG+CCAeuCklUSinVAK3etZpR748iLiKOty54i3h7fJOa+qM61SYMz1Baf/t3G2NWBPNmIjJcRNaKyHoRucvP8UkiskpEfhGRuSKS5XPsKhH5zbNoolJK1altB7Zx/rTzKXIU8c4F79AiqgURYRFNauqP6tTZnSUiYgOmYI24OhEY4xmi6+snoLcxpiswHXjMc25L4H6sEVp9sWbObR5tQKVUvdtfsp8L37+Qzfs389rI18hMzCQtLq3ZtCy86vJWxL7AemPMRmNMGTANON+3gDHma+/lL2Ah4H0W6DDgK2PMXs+U6l9h3QuilFIhVeIs4Q8f/YHF2xfzzPBn6JrSlZTYlEb9qNWjFVTCEJGhR9quRjqw1Wc7x7OvKuOA/x7luUopVWMOl4MbP7+RWetnce/Aezmzw5kkRSc1m07uysQYE3hhkWXGmJ5VbVdz7sXAMGPMdZ7tsUBfY8yf/ZS9ApgIDDLGlIrI7UCkMeYRz/F7gSJjzD8rnTceGA+Qmpraa9q0aQF/tsoKCgqIi4s76vNDReMKjsYVHI2rovd+f49XNr3ChWkXMr79eGxiq9Bn0RR/X0OGDFlqjOnt96AxJuAFWFZp+6cgzj0ZmO2zfTdwt59yZwKrgRSffWOAl322XwbGHOn9evXqZWri66+/rtH5oaJxBUfjCo7GZXG5Xea5H58zPIA5Z+o55rc9v5mNezcap8tZr3EFqiZxAUtMFfVqQJekRORNEXkDaCcib3jWg7UYOE5E2otIBHAZMLPS+/TwJIORxpg8n0OzgbNEJMnT2X2WZ59SStUqt3Hz0aqPmDR7Ev3S+/HPs/6JMYb0hPQmfyd3dQKdrfYtz8/TgLeP5o2MMU4RmYhV0duAN4wxK0XkIayMNhN4HIgDPvSMPvjdGDPSGLNXRB7GSjoADxlj9h5NHEopVRVjDPM3z2fczHG0T2rPq+e9CkBGQgYRtoh6jq7+BZQwjDHzAUTkoHfdeyiYNzPW7LazKu27z2f9zCOc+wZwNC0bpZQKyM+5PzPmozHEhsfyrwv+RbgtnLT4tGYx7UcggnoeBlBWzbZSSjVKm/Zt4uLpF1PkKOLjSz8mMSqR1rGtm/R05cEK9gFK/Y+0rZRSjVFeQR6XTr+ULfu38O6od2mX2I7EyMRmO3y2KsG2MJRSqknJL8nnqhlXsXj7YqaMmEKvtr2ItEfSOrZ1s7uTuzqaMJRSzVaRo4hbvriFLzZ8wV9O+wvDjx1OGGG0jWuaz+SuKf2NKKWapRJnCQ/Pf5i3fn6LcT3GMa7HONxutw6fPYJAHtHaMoDXcRtj9tdCPEopFXKlzlJeXvIyf//u74w4bgR/Oe0vlLpKyWqR1axmnw1WIJektnuWI13MswHtaiUipZQKIYfLwUerP+KOr+6gX3o/nhn+DMXOYjISMoiyR9V3eA1aIAljtTGmx5EKiMhPtRSPUkqFjNPtZO6muUz4bAJZLbJ4beRrON1O2sS10eGzAQgkYZxcS2WUUqreuNwulmxbwjWfXkNseCxTR00lPCzcmn02WofPBqLahGGMKamNMkopVV/cxs2a3Wu4asZVFJQV8PGlH5MUlUR0eDStY1rXd3iNRsCjpETkYhGJ96z/RUQ+FpGApjZXSqn6Yozh9/2/M27mODbt38RrI1+jQ4sO2MPszep53LUhmGG19xpjDorIAKwn4L0NvBiasJRSquaMMeQW5HLjrBv5cduPPDXsKfql99PZZ49SMAnD5fl5DvCiMeZTQKdvVEo1WHuL9zJ57mRmrZ/F5NMmM7LTSEqdpWQkZujw2aMQTMLYJiIvA5cAs0QkMsjzlVKqzuwr3scT3z/Bmz+/ybXdr+WPvf5IYVkh6fHpOnz2KAVT4V+C9SyL4Z6b9FoCt4ckKqWUqoEDJQd4++e3+cd3/2DEsSO4f9D9FDoKSY1LJS6y4T1StbGoNmGIyDIAY0yRMeZjY8xvnu0dxpgvfcsopVR9KywrZMbaGdw55076pPfh2bOfpdhZTHJ0sg6fraFA7sM4QUR+OcJxARJrKR6llDpqxY5ivt70NX/+75/JSszijZFv4DZu4iLiaBXTqr7Da/QCSRjHB1DGVX0RpZQKnVJnKUu2L2H8Z+OJtkfz7qh3ibJH6fDZWhTIjXtb6iIQpZQ6Wg6Xg9W7VnP9f67nYNlBPr70Y1JiU3C5XaTFp+nw2Vqio5yUUo2a0+1k/d713DDrBjbs28Cr571Kp+ROlDnLyEjQ4bO1SROGUqrRcrldbM3fyp1z7mRhzkKePOtJBrQbQFFZERmJGUTaI+s7xCZFE4ZSqlFyG+tcfiwAACAASURBVDfbD27nrwv+yn/W/Yd7BtzDqBNGUVBaQNv4tsSEx9R3iE1OMHNJiYhcISL3ebbbiUjf0IWmlFL+GWPYWbCTl5a8xOs/vc413a/hT33+REFZAa1iWpEYpQM3QyGYFsYLWNOYj/FsHwSm1HpESil1BMYY8grz+GDlBzz67aOcfezZPDj4QYocRSREJpAck1zfITZZgQyr9epnjOnpfViSMWafiOhcUkqpOrWnaA9fbviSO+fcSe+03jx39nOUucqIsEWQGpeqw2dDKJgWhkNEbIABEJHWgDuYNxOR4SKyVkTWi8hdfo4PFJFlIuIUkdGVjrlEZLlnmRnM+yqlmoZ9xfv4IecH/vzfP5OZmMmb579JmIQhCOkJ6YSJdsuGUjAtjGeBT4AUEfkrMBr4S6Ane5LNFGAokAMsFpGZxphVPsV+B64GbvPzEsXGmO5BxKuUakIOlBzgl52/MOHzCUTZo5g6airxkfGUOErIapGFPSyY6kwdjYB/w8aYqSKyFDgDazqQC4wxq4N4r77AemPMRgARmQacDxxKGMaYzZ5jQbVclFJNm9u4Wbt7LTd8fgMHSg/w0SUfkRafRkFZAe0S2+nw2Toixpi6eSPrEtNwY8x1nu2xWP0iE/2UfQv4zBgz3WefE1gOOIG/G2Nm+DlvPDAeIDU1tde0adOOOt6CggLi4hrerJYaV3A0ruA0xLjcxs2+A/t4ZNMjrDiwgkc6P0KvpF643C4ibBH1ehmqIf6+oGZxDRkyZKkxprffg8aYgBagN9YlqWXAL8AK4Jcgzr8YeM1neyzwXBVl3wJGV9qX5vnZAdgMHHOk9+vVq5epia+//rpG54eKxhUcjSs4DS2uAyUHzMq8lWbIlCGGBzDPLHzGbDuwzazZtcbsKthV3+E1uN+XV03iApaYKurVYC76TcV6/sUKguzs9sgBMn22M4DtgZ5sjNnu+blRROYBPYANRxGHUqoROFBygG0HtvHUwqf4etfX3D3gbkafOJqC0gIdPltPgkkYu4wxNRmdtBg4TkTaA9uAy4DLAzlRRJKAImNMqYi0Ak4FHqtBLEqpBiy/JJ+t+Vt5ZMEjvL/yfS5Iu4Ab+9xIsaOY6PBoHT5bT4JJGPeLyGvAXKDUu9MY83EgJxtjnCIyEeupfTbgDWPMShF5CKsJNFNE+mBd9koCzhORB40xnYETgJc9neFhWH0Yq6p4K6VUI7aveB+b92/mjjl3MGfjHG7tfyvD7cMpc5UhCG3j2+rw2XoSTMK4BuvZGOGUX5IyQEAJA8AYMwuYVWnffT7ri7EuVVU+73vgpCBiVUo1QnuL97Ju9zpu+uImlmxfwt/O+BtXdbuKXxf9itPt1OGz9SyY33w3Y4xW2kqpkNhTtMe6z+KzCWzO38xL577EuR3Pxel24jZuMhIyiLDp5BL1KZiEsVBETtRLQUqp2mSMYU/RHn7c9iPjPxvPgdIDvHvhu5za7lTKXGWUOa1pP6LDo+s71GYvmIQxALhaRDZi9WEIYIwxXUMSmVKqyTPGsLtoN3M3zuWGWTcQHhbOR5d8RJeULhQ7igHIapHFdgl4QKUKoWASxjA8SSJEsSilmhHjmXX20zWfcvPsm2kT24Z/X/RvslpkUVhWSIQtgvSEdO2zaECq/UuIyLfGmAHASiomC2/ySAhRbEqpJsoYQ25BLu/8/A6T/zeZE1ufyL8u/BetYlpxsPQgCZEJpMal6mioBqbahOFJFhhj4kMfjlKqqXMbNzsO7OCZRc/w+PePM6DdAF4f+Tox4TEcLD1Iq5hWJMck630WDVAwT9z7RyD7lFKqKm7jZtuBbdw7714e//5xRnYayTsXvEOUPYqCsgLS4tNoFdtKk0UDFUx7b6iffWfXViBKqabN5Xaxad8mJs6ayJvL3+Ta7tcyZYT10M4SRwntEtuREKVXuBuyQPowbgD+BHQQkV98DsUD34UqMKVU0+Fyu/htz2+M/2w8C35fwF0D7mJin4mUOEsQhOykbL3HohEIZPjBv4H/Ao8Cvk/JO2iM2RuSqJRSTYbT7eSX3F+45tNr+HXXrzwx9AnGnDSGgtICosOjaRvfVkdCNRKBdHrnA/nAmNCHo5RqSpxuJwu3LuTKGVey4+AOXh/5OkM7DOVAyQESoxJ1JFQjE3BaF5FI4CIg2/c8Y8xDtR+WUqqxc7gczNk4h2s+vYYSZwnTRk+jZ9ueHCw7SEpcCklRSdq53cgE0w78FKulsRSf2WqVUqqyMlcZn6z6hOs/u57YiFg+ufQTOiR1oNBRSFpcmnZuN1LBJIwMY8zwkEWilGoSylxlvLHsDW6efTOZCZm8d9F7JMckU+YqIysxS+eEasSCuXj4vYjobLVKqSqVOkt5/LvHufG/N9K5dWdmXDaDltEtAWtOKE0WjVuwkw9eo5MPKqX8KXYUc8/ce3j6x6c5Pft0Xj7vZdzGTYQtgrT4NGxhtvoOUdVQMAlDb9JTSvlVUFrAhM8nMHXFVEafOJrHz3ycUlcpLaJa0Dq2tY6EaiKCSRhXVbFfR0kp1YztK97H5R9dzhcbvuCG3jdw16l3UeQsIiVWR0I1NcEkjEKf9SjgXGB17YajlGpMdhzcwaj3R7Fw20LuHXgv43qMo9hZTHp8OvGROl9pUxNwwjDG/NN3W0SeAGbWekRKqUZhw54NnDftPH7b+xvPDn+Wczqeg8PlIKtFFlH2qPoOT4VATe7HjwE61FYgSqnG4+fcnznvvfPYXbSbt85/i37p/QgjjMwWmYTbwus7PBUiwdzpvYLyByjZgNZo/4VSzc43m79h1AejcBs3749+n47JHYkOj6ZNXBsdCdXEBZQwxOq1uhHY4tnlBHYaY5yhCkwp1fDMWDODP3z8B5Kikpg6aiqpcam0iGpBSmyKdm43AwElDGOMEZGnjDG9Qh2QUqpheuOnN5jw2QQ6JHXgXxf+i/iIeFJjU0mKTqrv0FQdCWZw9EIR6VOTNxOR4SKyVkTWi8hdfo4PFJFlIuIUkdGVjl0lIr95lqqG+CqlQuAf3/2D62ZeR8+2Pfng4g9IiEwgIyFDk0UzE0yn9xDgjyKyBWuIbVB3eouIDZiC9eS+HGCxiMw0xqzyKfY7cDVwW6VzWwL3A72x+lGWes7dF0T8Sqkgud1ubv/qdp5c+CTDjhnGP8/6J1H2KNIT0nUkVDNUl3d69wXWG2M2AojINOB84FDCMMZs9hxzVzp3GPCV94FNIvIVMBx4r4YxKaWq4HA5uPbTa3l3xbtc3uVy/jLwL0Tbo0lPSNeRUM1UMPdhbKm+1BGlA1t9tnOAfjU4N72G8SilqlBUVsRFH1zEFxu+4OZ+NzOh1wTiI+N1JFQzV5fPRfQ3hML42XfU54rIeGA8QGpqKvPmzQs4uMoKCgpqdH6oaFzB0biCU1BQwH+++g/3/HoPKw+sZOIxExluH07OihzsYXZ+47d6i6uh/r6aU1x1mTBygEyf7QxgexDnDq507rzKhYwxrwCvAPTu3dsMHjy4cpGAzZs3j5qcHyoaV3A0ruB8OPtD7lt7HxsKNjBlxBSGZA+hbXxbEqMS6zWuhvr7am5x1WgKSRFpE0TxxcBxItJeRCKAywh8apHZwFkikiQiScBZnn1KqVqyMm8lE3+ayNYDW3nrgrc4o/0ZtGvRrt6ThWo4ajrn8OuBFvTc5DcRq6JfDXxgjFkpIg+JyEgAEekjIjnAxcDLIrLSc+5e4GGspLMYeMjbAa6UqrkPV37IKW+cQpm7jKmjptI/vT9ZLbKICY+p79BUA1KjS1LGmHOCLD8LmFVp330+64uxLjf5O/cN4I2jCFMpVYUyVxmTvpjElCVTOCnlJP6v3f/RvU130hPSsYfV5RVr1RjovwilmqlN+zYx+sPRLNuxjKu7Xc2tJ9/KzpU7yUzM1AceKb+CmXxwkp/d+cBSY8zy2gtJKRVqM9bM4OoZV+MyLp4/+3nO7HAmbeLasDdsryYLVaVg/mX0BiZg3f+QjjV8dTDwqojcUfuhKaVqm8Pl4NYvbuXC9y8kLT6NTy75hBHHjSC7RbZ2bqtqBXNJKhnoaYwpABCR+4HpwEBgKfBY7YenlKotW/O3MvrD0Szatog/nPQHbj/5dtrEt9FnbquABZMw2gFlPtsOIMsYUywipbUbllKqNn2+7nPGfjKWUlcpTw97mrOPPZu28W31MaoqKMEkjH9jzVj7qWf7POA9EYnFZz4opVTD4XQ7uWfuPTz+/eMcn3w8Tw9/mk7JnUhLSCPCFlHf4alGJpi5pB4WkVnAAKypOiYYY5Z4Dv8hFMEppY7etgPbuGT6JXy/9XsuOfES7j7tbtLj00mOSdZLUOqoBDNK6lbgQ2PMMyGMRylVC75Y/wVXfHwFRY4inhj6BCM7jaRtXFviIuPqOzTViAXzNSMBmC0iC0TkRhFJDVVQSqmj43K7mDx3MiOmjqBldEumXzydiztfTHaLbE0WqsaCuST1IPCgiHQFLgXmi0iOMebMkEWnlApYbkEul02/jPlb5jPqhFFMHjCZdontSI5J1udtq1pxNHd65wG5wB4gpXbDUUodjbkb5zLmozEcLDvIo2c8yugTRpOWkKZzQalaFfAlKRG5QUTmAXOBVsD1gT6eVSkVGi63iwfmPcBZ755FfGQ8H4z+gCu6XkF2UrYmC1XrgmlhZAG3eKcBEZFTRWSKMebG0ISmlDqSvMI8xnw0hv9t+h8jO47k/kH3k52UTVJUkl6CUiERTB/GXSLSXUT+gdWHsQn4OGSRKaWqNH/zfC776DL2Fe/j4SEPc1nny0hPSCc6PLq+Q1NNWLUJQ0Q6Yj3saAxWv8X7gBhjhoQ4NqVUJW7j5tEFj3L/vPvJSMhg2kXT6JPeh9S4VJ2OXIVcIP/C1gALgPOMMevh0D0ZSqk6tLtoN1d8fAWzN8zm7GPP5uEhD3NMy2NIjEzUS1CqTgSSMC7CamF8LSJfANOw7vRWStWRb7d8y6UfXcquwl3cN/A+rux6JemJ6UTZo+o7NNWMVDtKyhjziTHmUuB4YB5wK5AqIi+KyFkhjk+pZs0Yw2PfPcbgtwdjExvvXfQeE3pPIDspW5OFqnPBdHoXAlOBqSLSEuu523cBX4YoNqWatb3Fe7nykyv5/LfPGdphKH87/W90TO5IQlRCfYemmqmj6iUzxuwFXvYsSqla9mPOj1z84cXkFuRy16l3Mb7neNIS0oi0R9Z3aKoZ02EVSjUgxhieXvg0d865k5TYFN658B2GZA/RhxypBkEThlINxP6S/Vw942o+XfspQ7KH8I8z/sHxrY/XhxypBkMThlINwJLtS7j4w4vJOZDDbSffxo19btSHHKkGRxOGUvXIGMPzi57ntq9uIzk6mbfOf4uzjjlLH3KkGiRNGErVk/ySfK6beR3TV09nYLuBPH7W43Ru3ZnYiNj6Dk0pv+r0K4yIDBeRtSKyXkTu8nM8UkTe9xz/UUSyPfuzRaRYRJZ7lpfqMm6lattvB3+j5ys9+WTNJ9zS7xbeufAderTpoclCNWh11sIQERswBRgK5ACLRWSmMWaVT7FxwD5jzLEichngnegQYIMxpntdxatUKBSUFfDykpe5+6e7aRnTkjfPf5Ozjz1bH3KkaoXbDU4nGBOa16/LS1J9gfXGmI0AIjINOB/wTRjnAw941qcDz4v+L1JNwC87f+HFxS/y7op3KSgroHdSb5678Dm6tumqz61QQXO5rMTgdEJZGRQXQ0kJOBwgYu0zxlqvTWJClYoqv5HIaGC4MeY6z/ZYoJ8xZqJPmV89ZXI82xuAfkAcsBJYBxwA/mKMWeDnPcYD4wFSU1N7TZs27ajjLSgoIC6u4T0DWeMKTn3GVeoqZf7u+czcPpOVB1YSLuEMbDWQc9qcwzH2Y0hIaHh3bOvfMTh1EZcx1uJ2l/+sTKRiciguLiAh4ejiGjJkyFJjTG9/x+qyheEv11XOVlWV2QG0M8bsEZFewAwR6WyMOVChoDGvAK8A9O7d2wwePPiog503bx41OT9UNK7g1Edca3ev5aWlL/H28rfZV7KP7MRsbjv5Ni4+8WKOaXkM8ZHxfL/ge/19BaEpx2VMeWvB4bBaCiUlUFpqJQdvIrDZrMVuh7Bqep9XrJjHoEGDa72FUZcJIwfI9NnOALZXUSZHROxAIrDXWM2gUgBjzFJPy6MjsCTkUSsVgDJXGTPWzODFJS8yb/M87GF2Ts8+nctOuozTs0+nVUwrosOjdahsM+btX/AmBu9lpMqXj7xJISam9i8p1VRdJozFwHEi0h7YhjVl+uWVyswErgJ+AEYD/zPGGBFpjZU4XCLSATgO2Fh3oSvl3+b9m3l5ycu8sfwN8grzSItP4899/8ylnS+lU3In4iPjdf6nZqa6/gWwfnoTQwO80lalOksYxhiniEwEZgM24A1jzEoReQhYYoyZCbwO/EtE1gN7sZIKwEDgIRFxAi5ggmcCRKXqnMvtYtZvs3hh8QvM3jAbEWFgu4E8OOhBhnYYSkpcCrERsdqaaOK8o5EKCsovI5WUWAnDmxjCwqzEEB4OUU1gNvo6vXHPGDMLmFVp330+6yVY06ZXPu8j4KOQB6jUEWw7sI3Xlr3Gaz+9Rs6BHFrHtOb6ntdz+UmXc0LrE0iMTNTWRBNTVf9CWZmVGMrKYPv28v6FqKjq+xfqgtMZmmtZeqe3UkfgNm7mbJzDC4tf4LN1n+EyLk7OOJnbTr6NEceOIDUuldiIWGxhtvoOVdVAIP0LIlYysNshOrp8uy4vKRUUwM6dkJdn/fQuvtt5eZCV1Y0RI2r//TVhKOXHrsJdvPHTG7y89GU27d9EUlQSY7uOZUyXMXRv252EyAR94l0jdKT+BShPDHXZv+C9rFU5AeTmWj991wsLDz8/KgpSUyElBY4/HgYPhri4nUCLWo9VE4ZSHsYYvv39W15Y/AIfr/mYMlcZvdr24o+9/sjITiNJi08jLiJOWxONgDcpOJ3W8NTiYuun01l3/QvGwIEDR24JeNeLiw8/PzraSgSpqdClS/l6Skr5emoqJCQcPppqxYodQKda/0yaMFSzt79kP28vf5sXl7zI2j1riY+I56ITLmJMlzH0SetDYlQiUfYonbqjgfHtX3A6rZZCcXF5/wKUtxZsNoiMtCrh2lBaClu3wtKlSaxadXirwLtdUnL4uTEx5ZV9t25WAmjT5vBEEBfXvIfVKtVgGGNYvH0xLyx+gQ9WfkCxs5guKV14aPBDjDp+FBmJGcRGxGIP0/8i9c3lspJDUVF5x3NxsbXuvbHNX/9CTRUWwubN1rJlS/n65s1WR7c1SUa3Q+Xj4spbAL16VUwAvuuNaRhtZfq/QTUrBWUFTP1lKi8ueZGfd/5MtD2aEceNYEyXMZyceTJJUUnamqhjLld534J35JF38SaFsjLIyanYvxBbw4l9jYF9+/wnhC1bYNeuiuVbtoSsLOjXz/qZnQ0Ox0/079+D1FSr5dBg+Js/pBZowlDNwi87f+GFxS/w7xX/5mDZQTomd2TyaZO5+MSLaZfYjvjIeG1NhIAx5QnBmxTKyqxLOmVl/mdWtdnK+xdqOhrJ7bYuD/lLCJs3W30Mvtq2tRLBmWeWJ4XsbGvd39RfK1fm07598HHVWEkJtm2/Y8/ZjG3rZmxbN2Hbuhn71s3Ytv1OZGoKDP+p1t9W/4eoJqvYUczs3Nnc9dpd/LjtRyJtkZx1zFmM6TyGgdkDSYpOItoera2JGnC7q08IlXkTgt0OERE1v3zkdMK2bf4TwpYtFfsRbDbIzLSSQI8e5QkhO9vaX1t9HDUlxUXYcrYcngy2bsa2/Xdsu3ZWKG9sdlwpbXClpFHavR95WW1pGYLpajVhqCbDGMOqXav4csOXfLnxSxZsWUCho5DsFtbkf5d2uZQOSR2Ij4gn3BZe3+E2Ct6E4L1c5HBUvFzkbSH4TnkRitFHZWVhrFvnPynk5FRMTFFRVosgKwsGDbJ+tm9v/UxPt+Kqb1JwsEJCOJQMtm2x9u/dXaG8sYfjSm2LKzWdkt4DcLXJwNU2A1daJs6sDrjTMq0PHh6OCbPx+28LOT4EX4Q0YahGbduBbczZOIcvN3zJ3E1z2VloffPKSsxixHEj6CW9uPzMy2kZ05KY8Jhm35rwnSq78uJtJezYUbH/wHseVLxcFBFx9AmhtBT27oU9e6yf3sW77bt/zx7YvXtghfMTEqxWwUknwXnnlV82ys62Opbr+25ryd+PLWeLdcno903lLYVtv1sJIX9fhfLuiEicKWk4W6dR3Od0HG0ycKRm4GiTiSOzPa42mdYwL7t3ulpb+R/Dl4EwYx0OxT91TRiqUTlQeoD5m+fz5cYv+WrDV6zdsxaAllEt6ZPeh1MzT+W0dqfRqVUn4iLiWPLDEjJbZFbzqo1PdRW/753L3sWbEI5UkTid1gikYEYbGQP5+eWV+759h1f8lfcXFPh/LRFISrI6mJOToUMH6N0bbLZN9O3b/lBSSEqq+yGnh/pjHG5k9y5kWw6tfvoG+9xlhOdsInzbZuzbtxC+fQu2goqdI+7IaByp6Thbt6VkwDCrddA2A3d6JiarPaSlERYdSViEnbAIO+HhNiLDbYgt7NAosEAXsFpdoaAJQzVoDpeDH7f9yFcbvuLLjV+yeNtiXMZFlD2Knm16Mqn/JAZmDaRn256H7r6OsEUcakmI30esNByhqvjDwsovDwXTX+Atv3t34N/+9+3z31cBVgvEW/m3bGldGkpKKt+u/LNFC+uLc2UrV26hc+fQ9i57kwF5eYRtz8G+Y6tn+R37jq2E79yGfWcO9rwdiKMMgCzPuSY2FtM2Ddq0xfTogis9HTIykMxM5JgOhLVtQ2RkpHU9zG4vvzmkvptCQdKEoRoUbz/EVxu/4quNX/HNlm8oKCsgTMI4sfWJXNP9Gga0G8ApmaccesZEhC2iXmeG9X0Smm/l7113OitW/N5houvXh67iB+uyz/795Ut+vlW5++6rvOzePYCiIv+vJ2JV6N4KPjvbut/gSAmgoQw1NcZKBiZ3J2Hbc7Bt34pt+++E5249PBk4HRXPDQ+H1FQkNRW6dYY2Z1jDqdq2ZYkx9B40CGndGmkCCaE6mjBUvdt+cHuFfojcglwA2iW2Y8SxIzg181QGZQ0iIzGDmPAYIu2RtZog/FXylb/1+1b23nXvncZVvWbljmDfG8xEArvkY4x1A9mePUeu6P0t/qab8AoLsyp/79KqFRxzDBiTy3HHZfhNAImJVl3Y0BgDbocLdu5EtuUQtn0r9u2/l7cQcrcRnrcNe952xFXxD2YiIg7dUSc9ToI2Qw8lA9LSoF07pE0bq/8gIqJiIrDZKPjmG6sjpZlogH9+1dQdLD3IvM3z+GrjV8zZOIfVu1cDkBSVRL+MfpyScQoDswZyYusTiY2IJdIWWe38TZUrfN8lP7/ipRxvxe9d9/Kt5Kur8L11RiCdvsZYlXd+Phw8aC2//NKSdevKK/d9+6zj/ir+qpISWPVYUpJV6ScmQrt20LVrxWTgb4mL8//ld+XK9XTunFH9hwox76U6Y8A4XUTs2oVZ+CNh27da9x/s2Io9N4fw3ByrZbAr9/BkEBkJKdZt1tLTJxmkpVlLZqbVavDXMmiCrYPaoAlDhZzD5WDRtkV8ueFL5mycw6Lti3C6nUTaIumV1ov/O/n/GNhuIL3T+hAXnoBdIrCHhVsVvgOKSytW8L7X8r3f+KFiJX/ovR3WnD7e//veSzsi1hfGyMgjf8t3uawO2oMHrZu8vBW+v23vvoKCw4955zYq17XCVnx8eaXfooU166i3cvcmBN/j3qU+7huoUJlXsVQu53suAG43tvy92HbvxL47F9uuHdh352LftQPbrlzse3Zi370T2948bPt2067SncsmMhJS20BKa6R3d2sypjZtypNBRkZ5MoiI0GRQSzRhqFrnchlW5q1mzsY5zNn0Jd/8Pp9CRwGCcHxyFy7rdA19UwbSK2UALSJaY5dwwkwEBblQgP+KHypez/edO8jft3xjrGv4+fnh5OUducL3reArV/j+ppOuzG63KvyEBOtnfDxkZJTvi4srP+bd3rt3GT169CQpydoX6nsD/HWu++7zlnG7qx7B5OWbdL3rYWEQJgZbQT72PTux7co9tITl7SBs5w4kbydheTuRvJ2wexfip+lkIiKgVSukZUtIS4aunaBVK9bZ7XTs3Nm6kSIjA0lJOTwZeJt+/nrNVa3QhNGEHekbXyD7vBVJ5Us8/vbtKNjNP+f8iwXbvuKH3DnklewAIC0mi4Gp59Iz+TT6pgwiPT6T6IhIIu0RuN1CSTEUlpRPKFdcbE0y5/3pu/ge991Xubx33Yr11CP+jqKjD6/U27Q5vIL3TQbexbsvKir4IZ4rVx7g2GODO8f3916538XfN3nfmHxnbfVeffEmXO+6iDWpXlZWpaGahQVI3k5kZy6yM9eaknXHDuvnzp0Vp2gtLT08aLu9vCMkORk6ZFudJq1aQevWFaZrleTk8iTgkwi2//ADHQcO1GRQzzRhBKlypVp539GU8S1buaKu/NPptP6v+qvAfSsRf3EH0sHqew3f98lj3mcK5B3IZ/OeHLbs3cH2/Dy2799D7v597DtYBI5YIlz9aWUbw4lhGcTTFpszkbwSGzOKhX8XSYVK3+E4cjyVhYdbo26io60lJsZa4uOtOsd7zLs/OhoKCtbRsWPHwyr6uDhrvT7u+vX+fbw3xvn+3SvPqwQV/y7eyt1745xvhe/9ku3v23+Fv31pqdWLvnO3NX52165DPzsuX07ks89UTAL+mllhYda1slatrCTQrVt5EkhJsRJBaqp1eSg5ubyfwLclUGWAfngznqpXmjD88E59sHHj4ZW2P8FWxkfi29FaubyIFUdRkRWf9/nCRUUVH0LvreR9v30HsxQVGYqLwRh/ASd6ls5+YjfYYwxl0UJBLJgYOVS5t0wq4sSNTwAAEEtJREFUr8j9Vfr+KvvK+46mcl+5cjudO3cM/kQflVthvsnZe7y6yUF9//5hYda27zd834q/8uJb+R/G7bZ6zHd6Kn3fJS+vPBl4b6zYvdu61laFlIQEq8JPTrY6UgYMsCp/7+JNAr79A77B+iYC1eRowvDDO/9+WJj1/9H7pC7fOXR8p2Cuar/DYZ1bef6d6socqVxx8Wl+W/3VCQs7vCIOj3RgwotwRR3AFbMPp+ymTHIpNrkYewGEF0F4ETExYbSKjyclMYG2LZLJTEqhfat0jmudRVJ8FLGxYeTmLKJ794H1NvVG5Rab77rDcfg3eN/1QJK53V5ecYeHV7xc7q34Kw+brWodrDtxMyvfgO59VmdepW/+lVoB7NpVXvnv21d1toqJqdhr3rmzte57K3VystUqSE2F1q35btUqBvftW7Hi964382lVlCaMw1gPUIeyskG1PqW83W5VNt7WeUTE4dsREdY8//72h4dDYeE2MjPbVaj4o6LKv7H7fnOPjgZ7pIPdjk1sKVrDxvy1rN+/hvX717Ahfx37SssnOAsPiyA9th3HxrenXfwxZMV3oH1CVzomdaZVZBtsYXZsYicMOyCHVby77e5qO4j9VcyB7quO7xda77pIeeL3vfHNd71yZV5VZR8w7xja/fnWONkqluNWrIAXX6yYEPbssb4d+GO3l1f+3lume/SoWPm3bFl+SahNG+u6W1XNlqoSwNq11rU6pfzQhFFJXBzccAPs2PE7GRlZFSrryhV4MPu9l29rauXKjXTu3O6w/f/f3vkHWVmdd/zz3XuXhd2FXViUAUR+BNSKOmgQaewQYm0aNSP9QUcyzVSntkxb07TTdlocE8f+YCZppjVtTWJJojExRoxNJoyTRtPiDq0KigEBBRQFC9UKsrusC7vLsvv0j3Pe3cvLvXfvst57Dft8Zt55z3vO++PLw93zvOfH+5y2nqO83rGX1zqCU3j94B72HdvDwffe4JQNzUZpHtfCrIYPcc3U65nVOI/ZjfOZM2kBcyddQn12AtmaLDXKko1i8/U05NZByVv2oUOhpyI90Jq7LzWv1PPT6Xy8+WaeN/lCDAyE/voiFX3RrbMzbMU+nIic39gYKvdkHGDBglDhJw4gcQLTpoXKf8qUoaZNvorf3/6dCuAOI0V9PaxdC5s27efyy2cPf0ER0rNXkq+D881IyndN+k37ZH8f/3P8IPv2bGD/e3vZ37mH1zv3cKDzVdpPnt5amNEwmwsb53Lt9OuY0zSPec3zWXjeQmZMms64bJbaTJZx2Sw1NTrjbTrdfVIKZ7vAzagxCwM1x4/n3c7fvBleeaX0yj7fqHMuNTVDo+fJyHlLS2iWJsfpaVTJxxOJM2hu5plt21i+dGl+B+A4H1DcYRQg6U4ezYB2uuUf0sZJO0HHyTY6+9o5drItbu109LbRGfcdvW209bTR0dNGR0877T3tdJ48dtr9Wya0MLd5HtfPv455zXOZP2U+F0+9mEumXjK4HnVoMVS5EkpG6gtU6nR1FS4rVH7iRMgvFPgocmnuQW3tmZX69Onh7T7XASRbU1Oo9JN+/8mTQ156OlLuqHSpb/vJoJLj/BxRUYch6RPAPwEZ4Btm9oVUeR3wbeDDwFHgFjM7EMvuBG4H+oHPmtmT5dKZTP6YOTPRdWb44H47RUdPOx297bT1tNE+WLGHfVt3G23dbbR3t3O0+yjtPe20d4f8voHC80mzNVma6prCNr6J5gmTmDP5QibVTaK5rplxHeP4yKKPsPD8hcyYOIPaTC3ZmiwZZYoPOCcR75I5sr29hdPFygqkLzt4MBguX4VfLKhRXiNkz5wmleyTz5uLTadKtsZGthw5wjVLloQ3+4aG0wdw8+0dxylIxRyGpAzwFeBXgEPAC5I2mNkrOafdDrSb2XxJq4AvArdIuhRYRZjLOQP4D0kXmdkZARdGy/GTx3lk5yNsPbSVpwZ+Eiv9NtpOHB2s8Nt62nmvr/jnsBMz9TRl6mnONNBcU8/8mgYmZ86jqXE8zYRtso1nykAdU6yOFpvA1IHxTOrPkukFnehHfX3oVD/q70d9p1B/P28d2MOMnxyG3vWhIi5UkSfrZCbHI/3ooRA1NUMeta5uMF2XLGRQXx+6aNKVd76KvaEhvNWn3+zHj8/fR18o+H8Rultbw+wgx3FGTSVbGEuAfWb2BoCkR4EVQK7DWAHcE9OPA/cpvDavAB41s15gv6R98X7Pvd8ie945xOonVgNQuw+mdA9tM7vh8m6Y3DOUNznZ9wwdN/dA7cAJ4ATwbtHnFSV3on4mrLTVknRlJBV27kh7Y+NQOrdST5+bU9EPLpuWm19fH/LS24QJoVsnTxfMi88+y/Jly3wA1nHOYSrpMGYCB3OODwHXFDrHzE5JOga0xPzNqWtnph8gaTWwGmDatGm0traOWKROHGfL7hsY35dhYnZCrKyzWCaDNWSxSZmQTm0DmQwnMxn+L5Ph7Wx2qCydrqkZzBvIZk8vT11zRheJRFdXF43p0eVyVNB9fWEr8pFXLl1dXbRu2vT+6xglXV1dZ/U7KDeua2S4rpFRLl2VdBj5arX0lJRC55RyLWa2DlgHsHjxYlu+fPkIJUZuvInW1lauONvry0hraytn/e8qI65rZLiukeG6Rka5dFVylO8QkDsj/gLgrULnSMoSYlC0lXit4ziOU0Yq6TBeABZImitpHGEQe0PqnA3ArTG9EthoZhbzV0mqkzQXWAA8XyHdjuM4DhXskopjEp8BniRMq33AzF6W9DfAVjPbAHwT+E4c1G4jOBXieY8RBshPAXeUY4aU4ziOU5iKfodhZj8GfpzKuzsn3QP8VoFr1wJryyrQcRzHKYh/qeQ4juOUhDsMx3EcpyTcYTiO4zgl4Q7DcRzHKQnZcOGcf06RdAR4cxS3mMqo4nqUDdc1MlzXyHBdI+Nc1DXbzM7LV3DOOozRImmrmS2uto40rmtkuK6R4bpGxljT5V1SjuM4Tkm4w3Acx3FKwh1GYdZVW0ABXNfIcF0jw3WNjDGly8cwHMdxnJLwFobjOI5TEu4wAEkHJO2UtF3S1pg3RdJPJb0W95MroOMBSYcl7crJy6tDgX+WtE/SDklXVVjXPZL+N9psu6Qbc8rujLr2SvrVMuqaJelpSbslvSzpT2J+VW1WRFdVbSZpvKTnJb0Udf11zJ8raUu01/oYTZoYHXp91LVF0pwK6/qWpP059loU8yv224/Py0jaJumJeFxVexXRVX57mdmY34ADwNRU3t8Da2J6DfDFCuhYBlwF7BpOB3Aj8O+ExaWWAlsqrOse4C/ynHsp8BJQB8wFXgcyZdI1HbgqpicCr8bnV9VmRXRV1Wbx390Y07XAlmiHx4BVMf9+4A9j+o+A+2N6FbC+TPYqpOtbwMo851fstx+f92fAI8AT8biq9iqiq+z28hZGYVYAD8X0Q8CvlfuBZraJENa9FB0rgG9bYDPQLGl6BXUVYnD9dTPbDyTrr5dD19tm9rOYfg/YTVi6t6o2K6KrEBWxWfx3d8XD2rgZcB3weMxP2yux4+PAL0vv/3rARXQVomK/fUkXADcB34jHosr2yqdrGN43e7nDCBjwlKQXFdYFB5hmZm9DqACA86ukrZCOfGukF6uUysFnYhP3AQ112VVFV2z+X0l4O/3A2CylC6pss9iNsR04DPyU0JrpMLNTeZ49qCuWHwNaKqHLzBJ7rY32uldSXVpXHs3vN18G/hIYiMctfADslUdXQlnt5Q4jcK2ZXQXcANwhaVm1BZVASeucl5GvAR8CFgFvA/8Q8yuuS1Ij8G/An5pZZ7FT8+SVTVseXVW3mZn1m9kiwjLHS4BfKPLsqumSdBlwJ3AJcDUwBfirSuqS9EngsJm9mJtd5NnV1AUVsJc7DMDM3or7w8APCX9I7yTNtrg/XCV5hXRUdZ1zM3sn/pEPAF9nqAulorok1RIq5e+a2Q9idtVtlk/XB8VmUUsH0Ero026WlCymlvvsQV2xvInSuyZHq+sTsWvPzKwXeJDK2+ta4GZJB4BHCV1RX6b69jpDl6SHK2GvMe8wJDVImpikgY8Duzh9ffFbgR9VR2FBHRuA34kzIJYCx5JumEqQ6gP9dYLNEl0VWX899g9/E9htZv+YU1RVmxXSVW2bSTpPUnNMTwCuJ4yvPA2sjKel7ZXYcSWw0eIoagV07clx+iKME+Taq+z/j2Z2p5ldYGZzCIPYG83st6myvQro+nRF7HW2o+XnygbMI8xQeQl4Gbgr5rcA/wm8FvdTKqDle4Suij7CW8HthXQQmplfIfRB7wQWV1jXd+Jzd8Qf5PSc8++KuvYCN5RR1y8RmtY7gO1xu7HaNiuiq6o2A64AtsXn7wLuzvkbeJ4w2P59oC7mj4/H+2L5vArr2hjttQt4mKGZVBX77edoXM7QbKSq2quIrrLby7/0dhzHcUpizHdJOY7jOKXhDsNxHMcpCXcYjuM4Tkm4w3Acx3FKwh2G4ziOUxLuMJwxhaS7FCKi7ogRPa8p8/NaJY1qbWVJN0taM8p7LE+imjrO2ZId/hTHOTeQ9IvAJwmRZHslTQXGVVnWsJjZBsJ3G45TVbyF4YwlpgPvWgidgJm9azEsjKS7Jb0gaZekdUmU0dhCuFfSJoX1La6W9AOFtRD+Lp4zR9IeSQ/FlsvjkurTD5f0cUnPSfqZpO/HWFPpcz4r6ZV4n0dj3m2S7ovp7Tlbt6SPxmgFD0T92yStKPDvnyTph/H+90vyv39nRPgPxhlLPAXMkvSqpK9K+mhO2X1mdrWZXQZMILREEk6a2TLC2gc/Au4ALgNuk5REI70YWGdmVwCdhLURBomtmc8B11sIdLmVsJ5BmjXAlfE+f5AuNLNFFoL0fT7e41nCV+Ibzexq4GPAl2KYmzRLgD8HLicEQfyNvFZynAK4w3DGDBbWXPgwsBo4AqyXdFss/pjCKmk7CUHmFuZcmnQH7QRethDkrRd4g6GgbgfN7JmYfpgQHiSXpYSFkp5RCON9KzA7j8wdwHclfRo4laccSQuALwG3mFkfIf7ZmnjfVkKIigvzXPq8mb1hZv2EcC9pjY5TFB/DcMYUsbJsBVqjc7g1dv18lRBj56CkewiVbkJv3A/kpJPj5G8oHWMnfSzCOg+fGkbiTYQVDm8GPi8p13ElATIfA34/6U6L9/5NM9s7zL2H0+g4RfEWhjNmkHRxfDtPWAS8yZBzeDeOK6w84+LhuTAOqgN8CvjvVPlm4FpJ86OWekkXpfTVALPM7GnC4jjNQHqc40HgQTP7r5y8J4E/zhl3ubKAxiUK61HXALfk0eg4RfEWhjOWaAT+JYbSPkWIKrrazDokfZ3Q5XQAeOEs7r2b0Fr5V0KU3K/lFprZkdj99T0NrYT2OcJ63wkZ4GFJTYRWw71RGwCSZhOc2UWSfjde83vA3xLWadgRncYBTh+DSXgO+AJhDGMTYe0XxykZj1brOKNEYRnWJ+KAueOcs3iXlOM4jlMS3sJwHMdxSsJbGI7jOE5JuMNwHMdxSsIdhuM4jlMS7jAcx3GcknCH4TiO45SEOwzHcRynJP4fkrE3Q7U/omkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "colors = ['red', 'green', 'blue']\n",
    "\n",
    "plt.figure(figsize=(6,4))\n",
    "for i, method in enumerate(methods):\n",
    "    plt.plot(N, T[:,i], label=method.title(), color=colors[i])\n",
    "    plt.fill_between(N, T[:,i]-3*S[:,i], T[:,i]+3*S[:,i], color=colors[i], alpha=0.1)\n",
    "    \n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xlabel('Sample size b')\n",
    "plt.ylabel(r'Avg. runtime [s] $\\pm$ 3$\\sigma$')\n",
    "plt.title(r'Data: Unit Gaussian in $\\mathbb{R}^{10}$', fontsize=10)\n",
    "plt.savefig('/tmp/normal10d_runtime.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run a simple sanity check that all VR PH computations return the same result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample data\n",
    "x = f(200)\n",
    "\n",
    "# Compute l1 distance matrix and get max. pairwise distance = threshold\n",
    "D_l1 = pdist(x, metric='cityblock')\n",
    "thr_l1 = np.max(D_l1.ravel())\n",
    "\n",
    "# Run Ripser\n",
    "dgm_ripser = ripser(x, maxdim=0, thresh=thr_l1, metric='manhattan')['dgms']\n",
    "\n",
    "# Run Ours\n",
    "X = torch.Tensor(x).to(device)\n",
    "l, _ = vr_persistence_l1(X.contiguous(),0, 0);\n",
    "\n",
    "# Run Dionysus\n",
    "filt = d.fill_rips(D_l1, 1, thr_l1)\n",
    "m = d.cohomology_persistence(filt)\n",
    "dgms = d.init_diagrams(m, filt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert(np.abs(dgm_ripser[0][:,1][:-1] - l[0].cpu().numpy()[:,1]).sum() < 1e-4)\n",
    "assert(np.abs(sorted([x.death for x in dgms[0]])[:-1] - l[0].cpu().numpy()[:,1]).sum() < 1e-4)\n",
    "assert(np.abs(dgm_ripser[0][:,1][:-1] - sorted([x.death for x in dgms[0]])[:-1]).sum() < 1e-4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}