Tuto-GUDHI-bayes-tda.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Bayesian Framework for Persistent Homology\n",
    "#### Authors: Christopher Oballe, Vasileios Maroulas\n",
    "\n",
    "Herein, we provide a tutorial for a Python module that implements the model for Bayesian inference with persistence diagrams introduced in https://epubs.siam.org/doi/pdf/10.1137/19M1268719.\n",
    "\n",
    "The module is hosted on github at https://github.com/maroulaslab/bayes_tda_python. It can also be installed from PyPI with the command:\n",
    "```shell\n",
    "pip install bayes_tda\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Usage\n",
    "\n",
    "We'll demonstrate how to use our code to create posterior distributions for persistence diagrams by looking at one of TDA's popular toy datasets, the noisy circle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "# noisy circle sample\n",
    "n_pts = 100\n",
    "noise_lev = 0.1\n",
    "t = np.linspace(0,1,n_pts)\n",
    "x,y = np.cos(2*np.pi*t), np.sin(2*np.pi*t)\n",
    "circ = np.array([[h,v] for h,v in zip(x,y)])\n",
    "noise = np.random.normal(scale = noise_lev,size = [n_pts,2])\n",
    "n_circ = circ+noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAD8CAYAAABpXiE9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXl0HNWV/z9XsiRL3t3CeIFuGTCLnRlCohD2sE0GnAyEH+QHiSCCAApmGOCXyRxMlJMEJoKQ5SSBkIBjHJZWyEKGsMcTYIhZQkBOMMZ4AJtYssGAbRbbMdiWdH9/VLXdknuprq7uqm7fzzl1urvqVdXtWr513333vRJVxTAMo1BqwjbAMIzKxMTDMAxfmHgYhuELEw/DMHxh4mEYhi9MPAzD8EUg4iEiC0TkLRF5IcvyY0XkPRF5zp2+HsR+DcMIjxEBbedW4MfA7TnKPK6qnw5of4ZhhEwgnoeqLgLeDmJbhmFUBkF5Hl44XESWAK8DX1HVZZkKiUgH0AEwatSojx544IFlNNEwdi8WL168XlX38LNuucTjL0BCVTeLyGzgd8CMTAVVdR4wD6C1tVV7enrKZKJh7H6ISK/fdcvS2qKqG1V1s/v9QaBORJrLsW/DMEpDWcRDRCaLiLjfD3X3u6Ec+zYMozQEUm0RkTuBY4FmEVkDfAOoA1DVm4AzgDki0g+8D5yl1p3XMCqaQMRDVT+XZ/mPcZpyDcOoEizD1DAMX5h4GIbhCxMPwzB8YeJhGIYvTDwMw/CFiYdhGL4w8TAMwxcmHoZh+MLEwzAMX5h4GIbhCxMPwzB8YeJhGIYvTDwMw/CFiYdhGL4w8TAMwxcmHoZh+MLEwygJ3d3dtLS0UFNTQ0tLC93d3WGbZARMOV+9YOwmdHd309HRwZYtWwDo7e2lo6MDgLa2tjBNMwLEPA8jcDo7O3cIR4otW7bQ2dkZkkVGKTDxMAKnr6+voPlGZWLiYQROPB4vaL5RmZh4GIHT1dVFU1PTkHlNTU10dXWFZJFRCkw8jKLI1KrS1tbGvHnzSCQSiAiJRIJ58+ZZsLTKsNYWwzf5WlVMLKob8zwM31iryu6NiYeRlfQqSXNzM83NzUOqJ5XQqmLJaiVEVSM7ffSjH1UjHJLJpDY1NSmQcWpqatJYLJZxWSKRCNt8Vc38H5qamjSZTIZtWmQAetTn/RmI5yEiC0TkLRF5IctyEZHrRWSFiDwvIh8JYr9GsKQ/pdvb23epkqSTWhblVhWrVpWWoKottwIn5Vh+MjDDnTqAnwa0XyMgUsHP3t5eVJWBgYG867z99tuRblWphGpVJRNIa4uqLhKRlhxFTgVud92kp0VkvIhMUdW1QezfKJ5MT+l87LP33rTNnk3bMcfAtm0wYoQzvfmm8zl2LNTVlcji/MTjcXp7ezPON4qnXE2104DVab/XuPNMPCLC8KfxCGB/4B+AfYG9gL3dz2nAWKChrw8mTsy94fHjYY89oLkZ9twT9t0XZsxwpv32g732gprSxO27urqGNCVDtKpVlU65xEMyzNOMBUU6cKo29oQoF6ocP3ky+61dyxHAPwIHAQ1pRdbjqP9q4Glg68iRXNbZCaNHw5gxjocxMAD9/c60fTu8+y6sXw/r1jmfL78MDz0EW7fu3PD48dDaCoceunOaMiWQv5WqPnV2dtLX10c8Hqerqysy1apKR5yaRAAbcqot96vqhzIsuxl4TFXvdH+/BBybr9rS2tqqPT09gdhnDGP1arj3XnjkEXj8cefmxnEFnwOWAi/V19OzbRuvAO8PW11EGBwc9Ly77u5uOjs7Wd3by8enTePqc87hxHgcnnsOnn0Wnn/eER+AWbPgn/8ZPvlJOOYYaGwM4A8bmRCRxara6mtlv800wyegBXghy7JPAQ/heCCHAc942aY11QZDMpnURCKhs0C/O26crp8+XRWcqaVFtb1ddcEC/d33v6+JeFxFRBOJxI71KLI51lOT6ZYtqk89pfrd76qeeKJqQ4Nj38iRqqeconrnnaqbNwd+bAr5D4lEYsixqQYooqk2KOG4E+ehtR0nnnE+cBFwkbtcgBuBlTgPtVYv291dxCPoCzN9e/tOmKCX1tTo4pRYgD5dU6N/OfNM1eXLPW2r2FwJXwL097+rPvSQ6qWXqk6d6tje1KR61lmq99+v2t/vef/FUs35IqGLR6mm3UE8gr4wk8mkNjU26lGgd4J+4ArGYtBLQCf79ByKETcRySgeIuJp3y3xuH4C9I7Ro/X90aN3ekzXXae6fn3JvYIgvK+oYuJRwQR6Yfb360XNzfonVzDWg/4Q9GCfN65qMF6R3/+YSVjHNjbqHy+9VPXYY1VB++vq9PbaWj2ghF5BMeIXdUw8KphALszt21Xnz1fdd19V0FdA54A2Zkkt9ypOQXlFfreTV3ReeEFvGz1aN4MOgHaDHlQCr8A8DxOPSFLUhTk4qPrb36oeeKAq6HP19Xo6aE0O0ShEAIK8afx4MF6EVUS0GfRa0E2uiPwSdD8IzA6LeZh4RJJsHdBisVjui/OJJ1Q//nFV0HenTtUz6+uzikV9fb3GYrGCqx5hu+texCu9TAy0yxWRD0D1iitUN27cUbYYEbDWFhOPSJJMJjP2UM14Ya9fr/rFLzqnbto01fnzdZ94PGf1xO+FHra77uVmz1Rmn8ZGXXH00c4x2nNP1QULVAcGQv8/UcTEowrIe2EPDuqTHR26vqZGt4P+dOxYvXP+fFUtnYdQanfdy9O8qDLPPKN6+OHOZX7ccTo9i8BWQ+DTLyYeVUBOAXjtNX3tH/5BFfRJ0A8Nu5FL+UQtlbtetjjC4KATTB47VjeL6MWgYp7HDkw8qoBsAnDBHnuoxmK6RUTnZLnwKzGgV/YqRF/fDgF+FHRqhRynUmPiUQUMF4BG0J+NGOGcokMO0QPzuNyVFtALJRg7OKh/uuAC3Syib4K2TZoU+eNUakw8qoSUAMRBn6+vd07Pf/yH6gcfVF2wL9T/s3y56qxZqiKq3/xmWVPdo4aJRzXx+OOqkyapjhmjet99O2Znqpr4bYKNAqFXtTZvVv3CF5xb4KSTVN97rzz7jRgmHtXCvHmqdXWqM2aovvjiLovTqyaxWEzr6uoqKs4xnNCrWoODqjfdpFpbq/qP/6i6enV59x8BTDwqncFB1auu2vkUfOedvKtUWzUmVBYuVB0zRv8+YYKePGVKRXpyfjHxqGQGBlT/7d+cU9He7vRT8UDY2Z/Vxv3XXqurRXQj6Ccq1JPzQzHiYS99CpPt26G9HW64Ab78ZViwwBk42AP2Jvpg+debbuJQVfqAB4Hjsdc05MPEIywGBhzhSCbhmmvge98raCBgexN9sPT19bEWOA5nxKr7gROx1zTkxK/LUo4pqtWWogN9AwOq55/vVFW+853w7NhN8HKchnewew70fTcXpJrBYh7lo9gmxuQdd+iCMWNUQX84bpwmk0kTgSLJdfy8nq/h5SaCPi+i20aOVO3pKfdfKhsmHmWkmFaOZDKpV9XVqYJ+Ly1Xw0uTqwlMZvKJQyHna/gx/u0NN6gmEk7ezYoV5f1jZcLEo4x4beXIdLNf3NysCnprllTzbBd36AlVESafOBTdKrV8uerEiar77af65pul+yMhYeJRRrw8yebMmbPLRXt0Q4NuAV0EWu9BPNIvbsvpyE4+cQjk2D31lGpjo+qRR6pu21aaPxISJh5lJJ8XkEwmd7mgp4G+DroStNmDcAy/uC2nIzv5xCEwr+0Xv3Bul0svDf5PhIiJR5nJFX8YfjHXgj4OuhF0pnvhpi/3EvMwzyM7XkcbCyRedPnlzi1TRdVFE48IMdxL+E9QBT2LnWNvDL+Q813cFvPITdmCydu2qR5zjFOFef75qghim3hEiHQv4Xic0bznu1WMYi6uarhQq4I33lCdPFnfjsd1fGNjxQu6iUeESHkJMTfO8SLoKNA5c+aEbZoRFPffrwp6TZ5YSyWIvYlHxEgmk3rXqFG6DfSkKVMie+EY/pkH2g96RIYgdiVVM4sRD3HWjyatra3a09MTthmF8/vfw8knw9e+Bv/5n2FbY5SAWfE4961ejQIfAj5w5ycSCQB6e3t3WSeRSLBq1apymegJEVmsqq1+1g2kY5yInCQiL4nIChGZm2H5uSKyTkSec6cLgthvJNm8GS66CA480BEPoyr56rXX8q8NDewLpC74VMfEbJ3pqq2TXdHiISK1wI3AycBM4HMiMjND0V+p6ofdaX6x+40s3/wm9PbCz34GDQ1hW2OUiLa2Ns6+5RbuaWriCuDoqVOZN28ebW1tu81wCUF4HocCK1T1VVXdBvwSODWA7VYeK1bA9dfD+efDUUeFbY1RYtra2jj1lVcYOXo0iw45hLa2NiD7cAmzZ8+mpaWFmpoaWlpa6O7uDsPswAhCPKYBq9N+r3HnDed0EXleRO4Skb2zbUxEOkSkR0R61q1bF4B5ZWTuXKivtzjH7sTUqSw+5RR44AE+LUJLSwsA7e3t1KSNz9Lf388tt9xCb28vqkpvby8dHR2VLSB+I62pCfgsMD/t9znADcPKxIAG9/tFwKNeth1Wa4uvZrYnnnAar666qvQGGpEhmUzq2MZGfQX0rzgv5aqvr9fa2tqCuyGEAWE21QKHAwvTfl8JXJmjfC3wnpdthyEevprZBgedd6JOneoM6W/sNqSSAtvcTOIzPPZdSk1h908KWzxGAK8C04F6YAkwa1iZKWnfTwOe9rLtMMTDVz+Shx92DuVNN5XNTiMapLoj1IAuc5MCawoQj0r2PIqOeahqP3AJsBBYDvxaVZeJyNUicopb7FIRWSYiS4BLgXOL3W+pyNac1tvbmz3Qdc01MGUKnHtu6Q00IkWqBWUQ+DpwEPA5j+tW/JizflWnHFOUPI/0aUg15k9/cryO739/l21VSoqy4Z/0aq6ALgX9q4jW1tTkvIZisVgkrgcsPT04MsU8Mk073M1TT3VGmtq0Ke92opqibBRH+kPiiokTVUH/u7NTY7GYtwdQiJh4BEz6xZDt5IuI6qpVzsuSOzt32Ua2CyfsOq5RYrZsUY3FVD/zmR2zojweSzHiYe9tyUBbWxurVq1icHBwR1+F4cTjcbjlFufHhRcOWdbd3c2GDRsyrldtKcrGMBob4UtfgnvugVdfBbKf80q/Fkw88pAtW/Caq692xOPkk2GYwOR6y1i1pSgbGbj4YhCBn/8cqN63+5l45KGtrY158+aRSCQQEWKxGI2Njfy6vR1ef50/HnDALuvkeqJUdHTd8Ma0aXDCCc7bAAcHq/ftfn7rO+WYojaeR3oQ9C6cwX7GNDbuEvjKVseNxWIhWW6Undtvd0KKixapanRb3rCAaXlIicIo0C2gN2QJfFlLi6GbNun2hgbtHj06r2CEKSwmHiUmdXJTQnCGm4r8iRwpxlF90hjlIZlM6i9qa3UDzgj62R4gYT9oTDxKSKaT+0vQN9PSkKPQ5GZEi0Qioae7D5mjczTPht2MW4x4WMA0D52dnWzZsmXH73rgU8DdOCnJVRH4MgKnr6+P/wa241wv6fOHl8u2ftQx8cjD8JN4GDAaeBBnTMrU6FGGkU48HmcTsIih4jG8ebaSm3FNPPIw/CQeDwwAfwzFGqNSSDXPPoAzQHKCzF5qRTfj+q3vlGOKYszjcdA/WyuK4YFkMqn/PGWKKujlOTrCVWpri716wQPd3d10dnayvreXd4DvAV9NWx7FIfWNiDAwABMmwNlnw09+ErY1uxD6qxeqnVRfl1agDnh82PJKCG4ZIVFbC4cdBk89tcui7u7uih4Q2cSjAE6YMAGAxcPmV0JwywiRww+HpUth48Yds7q7u+no6KjoAZFNPAqgbdYs1orwVtq8igluGeFx2GEwOAh//euOWcNTAAC2bNmSs1Nl1DDxKIB93n2XwYMP3tFJzppqDU/MmuV8Ll++Y1Yh+R1Rrd6MCNuAimHrVli+nGlz57LqW98K2xqjkth7bxg1aoh4xOPxjO+zHV4FTlVvUl5KqnoDhP7QMs/DK319TuR8//3DtsSoNEScdxeniYfX/I4oV29MPLzyt785n+4bwQyjIGbOHCIew8eJyVYFjnL6uomHV1LiMX16uHYYlcn06fDaa9Dfv2NW+nCXq1atylgNiXL6uomHR5bdfz/bgBHxeKSCVkaFMHmy08e2wPcvRzl93cTDA93d3Sx56CHW4vRrqcQ2eSNkJk92Pt94o6DVvFZvwsDEg/xNYZ2dnYweGODttHlRCVoZFUJKPNau9bxK6ro855xzALjjjjuyVm9CwW+nmHJM5egY52UkJxHRP4I+mundLYbhhZUrVUH15z/3VDyZTGpdXd2Q662uri7wTnPYYED+8dIUFo/HGQ+8O2zdKAStjAohFbd4//0ds3J5vJdddhnbt28fsont27dz2WWXlcVcLwQiHiJykoi8JCIrRGRuhuUNIvIrd/mfRaQliP0GQaZEHRjaFNbV1cUEEd5LWx6VoJVRITQ2Op8ffADk79uS7aVh2eaHQdHiISK1wI3AycBM4HMiMnNYsfOBd1R1P+AHwHXF7jcIuru7EZGMy9K9ira2NmLjxlE/enTkglZGhTBypPPpikeUk78847e+k5qAw4GFab+vBK4cVmYhcLj7fQSwHpyxRHJNpY55ZBt8VkR2rVtOmaJ64YUltceoYgYHnZjHfvupqmZ9D3IqjpbtXcdBv/uHkGMe04DVab/XuPMyllHVfuA9IJZpYyLSISI9ItKzrsA28ULJlqWnquZVGKXhLadPdr7krx/96EfU19cPWVZfX8+PfvSj0tpXAEGIRya/f/jwZF7KODNV56lqq6q27rHHHkUbl4tsJzDby60Nwzep4OcVVwD5k7/a2tpYsGDBkPyOBQsWROqhFoR4rAH2Tvu9F/B6tjIiMgIYB0PSJkKhoOy9mpohqcWGURBurCMV+/CS/OUlfT1MghCPZ4EZIjJdROqBs4B7h5W5F2h3v58BPOrWt0KloOy9ceOGjARlGAUxTDwg+uKQj6LH81DVfhG5BCcoWgssUNVlInI1TjDmXuAW4A4RWYHjcZxV7H6Doq2tzdtJGz8e3nmn9AYZ1UmqZSVNPCqdQAYDUtUHcd6DlD7v62nfPwA+G8S+QmPChIJSiw1jCKngf3NzuHYEyG6fYeoZ8zyMYkh1iJsyJVw7AsTEwyt77ul4Hj5DNVEdh9IoEynxmDy5aq4FG8PUK9OnO0GvN9/c2UPSI1Eeh9IoE26V985HH6Xj4our4lowz8MrqRHECngzXOoJc/bZZ1d+KrKRFU+eRF8fTJrEld/4RtVcCyYeXkmNXZoajjAP6R2fshGFcSiN4vD88qbly+HAAyM9JmmhmHh4JeV5rFjhqXimjk/DsS79lY+nDm6qjngcdFCkxyQtFBMPD3R3d9MycyavAA9++9ueAlz5niTWpb868ORJvPWW01J30EGRHpO0YPz2qCvHVI6RxPKRPtLYr0BXZhhpLBPZeuwCmkgkAh8RyigvyWQy7znewaOPOj1qFy4csq6IhH4tUESv2tAFItcUBfFIv0CucBxQHT/84siAl+ENjcok07nNeZ6//W3nVlu3Ljyjs2DiUULSx134pCsex3ocv3T4E2bOnDmReeIY/inYqzzlFNX99w/H2DyYeJSQ9AtlAugA6Nc8eB7DMU+kesg3kM8QBgdVm5tV29vLbqcXihEPC5jmIT3A9Q7wHHBiTU3BAa6qGHbOAAp8i9vKlbB+PRxxRImtKj8mHh5oTA1eC/xp5EiOqqmh7bTTCtpGNbXv7+4U1GLy6KPO51FHlcGy8mLikYNUAlD6iNUPq1Lb3w9PPVXQtqqpfX93p5BxYFbffDNramupmTWrovuxZMRvfaccU9gxj0yBsVGgW0H13/+9oG1ZzGP34xcLFuhm0B9H+JxjAdPSkC0w9gCoTp/uBMMKIErt+0bpOWfSJFXQk3PlgIRMMeIhzvrRpLW1VXt6ekLbf0tLS8a+KV+Jxfjuhg2weDF85CMhWGZUAj8R4Vyc1wR8kDZfRBgcHAzHqGGIyGJVbfWzrsU8cpAtMHbot74FtbVw110hWWZEnm3bOKumhvsYKhxQPXEuE48cZAuMffaii+D449m4YAEtiUTFD+piFIanLvi//z0TBwf5dUPDkNkV248lE37rO+WYwo555OLJjg5V0GMiHAwzgsdz4PuMM1T32EO7b7010nEuLOZRfg6Mx3l69WoeAM5Om59IJFhVwIBBRmWRLQ425Ly/844zVmlHB1x/fXkNLBCLeYTAy2vWcDvOS2gmps23pK/qxlOy34IFsHUrnHdemawKBxMPn8TjcX4GNLDzbVap+Ub1kjfZr78fbrgBjjkGDjmkjJaVHxMPn3R1dfFqUxOPA5fhjCRdVcEwIyN5U9PvuQd6e+Hyy0Owrsz4DZaUY4pywFTVCZ61u4lAX47FIhcMM0pDzmS/o45SbWlR7e8Pz8ACwDJMQ2RwUPXgg1UPOCDnBWPZpbsBjz3m3FI/+EHYlnimGPGwakuxiMBXvwovvQR3352xiOcRto3KRRW+9jWnleVLXwrbmrJQVFOtiEwEfgW0AKuA/6uqu7yTUUQGgKXuzz5VPcXL9qPcVDuEgQGYNcvJOl2yBEYMfZeWp+Y9o7JZuBBOOgluvBEuvjhsazwTZlPtXOARVZ0BPOL+zsT7qvphd/IkHBVFbS1cey28+CLccssui20sjyon5XUkEnD++WFbUzaKFY9Tgdvc77cBnylye5XLZz4DRx8NX/86bNo0ZJGN5VHl/OIX0NPjnPth6ejVTLHisaeqrgVwPydlKTdSRHpE5GkRySkwItLhlu1Zt25dkeaVERH43vecd3Rcd92QRbNnz864Sm9vr/WJqXQ2boSvfAU+9jFob89fvprIF1EFHgZeyDCdCrw7rOw7WbYx1f3cByc2sq+XaG5FtLYM5/OfV62vV12+fMesXKNtY31iKpsvf1lVRPWZZ8K2xBeE1VQLvARMcb9PAV7ysM6twBletl+R4vHGG6oTJjjt/QMDqpp9UKH0qba21gSk0li6VLW2VvXCC8O2xDfFiEex1ZZ72Zmd3Q7cM7yAiEwQkQb3ezNwJPBikfuNLnvuCd//PjzxBMybB3iLbQwMDFjzbSWxfbvTd2X8eLjmmrCtCQe/quOIFjGcVpZX3M+J7vxWYL77/QicZtol7uf5XrdfkZ6HqpM4dsIJqmPGqPb15X3DWPoUpSHqjBzJfVdd5Tjuv/lN4etGCCzDNIKsXKk6apTqMceo9vcPebdprmqMlzfRGeUh29gdD159teqIEaptbQWvGzUBMfGIKrfd5hzib35zyOxkMqm1tbXmeYSIF68gU6C7CfTlujrVadNU33476/azBcmjdn5NPKLMOeeo1tQ4/R7SqJQnUzXi9dhn8hCT7itH9Q9/yLmPgl5JGSImHlFm0ybVGTNUp05VXbt2yKJKqBNXI169guHlLnVySfU748cHto+wMfGIOs89p9rUpHrooapbtoRtzW6PV68g3UM5GnQb6L21tZq8/fa8+6gUz9LEoxK4+24nmeisswp+WZQRLIV4BclkUo+eOlXfAF0xYoT+vy9+0bO3WAmepYlHpfDtb2umAKpRXgryCt58U3W//VQnTtR7v/OdivAmCsHEo1IYHFRtb3cO+09/GrY1uzWevIJNm1Q/9jHVkSNVn3yyYuIYhWDiUUls3ar66U87VRgPdWejfKQLyox4XNccfLDTUnbPPaqqOZP7KpVixMNGEisDqTeMiQgjmpoYef/9PNHQwGB7O/z2t2GbZzB0tLd6VX7Q18e0JUt4+rzz4BRnCJra2tqM62abX/X4VZ1yTNXgeWRLTW8CfbKmRvtra3OmOO/ulCvomKqSjAR90G2SvTCtWpJMJs3zsGpLecnVHX8s6DMNDU4VZt68nNuphMh90JSzuVNEtBF0oZsE9sXhYt/UpLFYzGIeJh7lI193/CZQnT3bORXXXpuxGbdScgaCppwByo/utZf+2RWO9iznKhaLVd15MPGIMPkGAkokEqrbtjmdrED/95/+SfeNx4d4GNUY5fdC2VK8ly3TTXvsoZtB/yXHuRKRqvMATTwiTK7u+EOeWgMDusz1QB4GnZhWJtfFXM3kEt7AbtxHHlEdN051zz312tNPz+spVoNgpGPiEXHSvYdUb9pMF2EikdAvgL4PuhL0Q+4Fu7v2wM03DkpRVYaBAdVrrnGaYmfOVF21Kq+XWC1VlXRMPKqE1FPvUNDXQDeDnpfFA0mVjdqTMGi3Ple1zbeArl+/M8505pmqGzeqqrfhIqtNuE08qoT0m2QK6KOgCnp/U5P++qabsg4mFJUnYSkDu4HFPxYtUo3HnUGqb7xxSIDaq+dRTVVGE48qYfjNVwPaWVenA7W1zuAzjzwS6eBpKW0rett//7vqZZc5zeLTp6s+++wuRTKJXzbRisLxDgITjypiuNs/Z84c/dTkyfq/rheyADQW0SdhKVtHivJqFi1S3Xdf53K/5BKnz0qO/Qw//tXWPJuOiUeVkn7DNIJegzOmxDrQc0v4JPQbtyi1V1SwXWvXqp53nuNt7LOP6v/8T3n2W0GYeFQpmW7GWaCPu17IE6BHBPwkLOYJn2nd+vp6jcVi5b3xPvhA/3LmmbpJRLeC/nTsWL1z/vzS77cCMfGoUrJWA0D/Y+JEfd0VkdWHHOK8gCgAik3BTn9Kx2IxraurK5/Lv3276m236cZJk1RB7wbdtwqrGkFi4lGl5K0GbN6s2tWlOnasqoi+euSR+skpU3w/5XN1/vITtyhbcHf7dtWf/9wZtAd0WV2dnljFQc4gMfGoUjxXIdav12WzZ+tm1xN5EPQ40KbGRk0mk57r7EHnU5Q8vfzdd1Wvv35nMPSQQ1TvvltrAhTAasfEo4op5MafCNoJ+oYrIn8B/cqoUTq1sdFT1SFXkpQfl79knsfSpaoXXeS8VAtUDztM9d57d+RsRLk5O2qYeBhDbvwG0AtAn3NF5H3QbtATQWtz3EjZbrpYLObLpkCTxl57TfWHP1T9+Medy7ahwWlJ6ekp7X6rHBMPI+uNfwjoDaAbXCFZh5Mrcio4iVNplOKmK6qZ829/U/3JT1SiBzOdAAAIU0lEQVSPO85pbgXVD39Y9bvfdVLMS7Xf3YjQxAP4LLAMGARac5Q7CXgJWAHM9bp9Ew/vZLvxU60nDaCng94B+o4rJNrYqHryyc44Ik89pbp1a9luuoz72bBB9b77nESuGTMcG0F1//1Vv/EN1eXLS2LL7kyY4nEQcADwWDbxAGqBlcA+QD2wBJjpZfsmHoWR6YbMJCpjGxv14blznZv0oIN23qSNjarHHqt6+eWqCxY4KdzDvJNA7LzjDt135Ej9BOjlbpVqRcqzAOcFWbNnO9WU5ctzvufGPIziKEY8RlAEqrocQERyFTsUWKGqr7plfwmcCrxYzL6NXWlra6OtrS3jss7OTvr6+ojH43R1dXFCerm33oInnoBFi+DJJ+Hmm+H9951lIrD33junvfZypnHjYPRoZxozBurrob/fmbZvd6Z334X1651p3TpYuxZWrOC0pUtpcx4sAKwGnlHlN+PHM/e//guOOAIaGvL+39SgxVu2bAGgt7eXjo6OHcfCKC2iaSfR90ZEHgO+oqo9GZadAZykqhe4v88BPq6ql+Tbbmtrq/b07LJJo9QMDMDKlbB0qTO9+iqsXg1r1jjTBx8Uvs0JE2DSJNhvP374wAO8jFOHXQq84RYREQYHBz1vsqWlhd7e3l3mJxIJVq1aVbiNuyEislhVW/2sm9fzEJGHgckZFnWq6j0e9pHJLcmqWCLSAXQAxONxD5s3CqG7u3sXL2SXp3RtLey/P93PPkvnrbcOLfv5z8Pbb8OmTbB5887PrVuhrg5GjNg5jR8Pzc0Qizm/Xb7V3MyGDRt2sa3Q893X11fQfCNg/NZ30idyxzwOBxam/b4SuNLLdi3mESyFtKaUqrkzmUzukrIOTh+YQrdt+RzFQ9hNtXnEYwTwKjCdnQHTWV62a+IRLIXcbH5vzHwBzGJySXa37vLlIDTxAE4D1gBbgTdTHgYwFXgwrdxs4GWcVpdOr9s38QiWQtLF/aSWe/FW/KasZ9v2nDlzrLWlCEL3PEo1mXgES6k9Dy/r+PVorIpSGooRD3tX7W5EV1cXTU1NQ+Y1NTXR1dVVVNkUXgKYfrbrddtGmfGrOuWYzPMInkKSqgpNwPLqHfhJ7DLPozRg1RYjCpSyQ5p1disNJh5GZChlurilogdPMeIRSIZpqbAMU8MoLcVkmFrA1DAMX5h4GIbhCxMPwzB8YeJhGIYvTDwMw/CFiYdhGL4w8TAMwxcmHoZh+MLEwzAMX5h4GIbhCxMPwzB8YeJhGIYvTDwMw/CFiYdhGL4w8TAMwxcmHoZh+MLEwzAMX5h4GIbhCxMPwzB8YeJhGIYvTDwMw/CFiYdhGL4w8TAMwxdFiYeIfFZElonIoIhkffeDiKwSkaUi8pyI2ItYDKMKGFHk+i8A/we42UPZ41R1fZH7MwwjIhQlHqq6HEBEgrHGMIyKoVjPwysK/LeIKHCzqs7LVlBEOoAO9+dWEXmhHAZ6pBmIkvdk9uQnajZFzZ4D/K6YVzxE5GFgcoZFnap6j8f9HKmqr4vIJOAPIvK/qrooU0FXWOa5++7x+x7NUmD25CZq9kD0bIqiPX7XzSseqnqi342nbeN19/MtEbkbOBTIKB6GYVQGJW+qFZFRIjIm9R34JE6g1TCMCqbYptrTRGQNcDjwgIgsdOdPFZEH3WJ7Ak+IyBLgGeABVf29x11kjY2EhNmTm6jZA9GzqWrsEVUN0hDDMHYTLMPUMAxfmHgYhuGLyIhHFFPdC7DpJBF5SURWiMjcEtozUUT+ICKvuJ8TspQbcI/PcyJybwnsyPl/RaRBRH7lLv+ziLQEbUOB9pwrIuvSjskFJbZngYi8lS1HSRyud+19XkQ+ErI9x4rIe2nH5+ueNqyqkZiAg3ASVh4DWnOUWwU0R8UmoBZYCewD1ANLgJklsuc7wFz3+1zguizlNpfwmOT9v8DFwE3u97OAX4Vsz7nAj8txzbj7Owb4CPBCluWzgYcAAQ4D/hyyPccC9xe63ch4Hqq6XFVfCtuOdDzadCiwQlVfVdVtwC+BU0tk0qnAbe7324DPlGg/ufDyf9PtvAs4QUrXh6Gcx98T6iRAvp2jyKnA7erwNDBeRKaEaI8vIiMeBZBKdV/sprKHzTRgddrvNe68UrCnqq4FcD8nZSk3UkR6RORpEQlaYLz83x1lVLUfeA+IBWxHIfYAnO5WEe4Skb1LZItXynnNeOVwEVkiIg+JyCwvK5SrbwtQ/lT3MtmU6Ynqu/07lz0FbCbuHqN9gEdFZKmqrvRr0zC8/N9Aj0kevOzrPuBOVd0qIhfheEXHl8geL5Tz+HjhL0BCVTeLyGzgd8CMfCuVVTw0gqnuAdi0Bkh/ku0FvF4Ke0TkTRGZoqprXTf3rSzbSB2jV0XkMeAQnLhAEHj5v6kya0RkBDCOErjNXu1R1Q1pP38GXFciW7wS6DVTLKq6Me37gyLyExFp1jxDaFRUtSWiqe7PAjNEZLqI1OMECANv4XC5F2h3v7cDu3hGIjJBRBrc783AkcCLAdrg5f+m23kG8Ki6kbkSkNeeYfGEU4DlJbLFK/cCX3BbXQ4D3ktVR8NARCanYlIiciiOLmzIvRaRam05DUeRtwJvAgvd+VOBB93v++BE05cAy3CqFqHapDuj5y/jPN1LZhNO3OAR4BX3c6I7vxWY734/AljqHqOlwPklsGOX/wtcDZzifh8J/AZYgdMlYZ8Sn6d89lzrXi9LgP8BDiyxPXcCa4Ht7vVzPnARcJG7XIAbXXuXkqN1sUz2XJJ2fJ4GjvCyXUtPNwzDFxVVbTEMIzqYeBiG4QsTD8MwfGHiYRiGL0w8DMPwhYmHYRi+MPEwDMMX/x+DNn/oLGsJYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(circ[:,0],circ[:,1], color = 'red') # circle\n",
    "plt.scatter(n_circ[:,0],n_circ[:,1], color = 'black') # noisy sample\n",
    "plt.xlim([-1.5,1.5])\n",
    "plt.ylim([-1.5,1.5])\n",
    "plt.axes().set_aspect('equal')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll compute the Rips filtration for our point cloud."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gudhi as gd\n",
    "\n",
    "st = gd.RipsComplex(points=n_circ, max_edge_length=3.).create_simplex_tree(max_dimension=hom_dim + 1)\n",
    "dgm = st.persistence()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this demonstration, we'll focus on modeling the distribution of the feature that corresponds to the large hole in the middle of the circle. Hence, we only consider 1-dimensional homological features. Let's compute and visulaize the 1-dimensional persistence diagram for our data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VdW5x/HvC0ghCAiCw2UKVaKgRaBRrBoFp4LXinORalFRnLD1autQvah4qaJWFGcUiloncCpFkBaLmKooYRAFSoqoGKkCihNzwnv/WCE9hgwnkH32Ocnv8zw85uyzc87PPEnerLX2fpe5OyIiIgAN4g4gIiLpQ0VBRETKqCiIiEgZFQURESmjoiAiImVUFEREpIyKgoiIlFFREBGRMioKIiJSplHcAWqqTZs2np2dHXcMEZH0t3UrfPoprFnD3K1b17h72+o+JeOKQnZ2NgUFBXHHEBFJby+/DJdcAqtXw7Bh2L33fpzMp2n6SESkLnGHwYPhxBOheXN44w0YMybpT1dREBGpC7Y1NzWDLl3g5pth/nz4yU9q9DIqCiIime7DD+GnP4Vp08LjG26A4cOhceMav5SKgohIpiouhrvuggMPhLfegi+/3OmXzLiFZhERAd59Fy64AAoKwvrBAw9Ahw47/bIaKaSRwsJCxo0bR2FhYdxRRCTdFRTAihXw7LMweXKtFATQSCFtFBYW0qtXL9wdM2PevHnk5OTEHUtE0smsWfD553DmmXD++XDaabDbbrX6FhoppIn8/HzcnfXr1+Pu5Ofnxx1JRNLFV1/BRRdBnz5w663hpjSzWi8IoKKQNvLy8jAzsrKyMDPy8vLijiQi6eDFF6FbN3j0UfjNb+Af/4AG0f3qjuyVzWy8ma0ys/erOe9gMysxs9OjypIJcnJymDdvHmPGjNHUkYgECxbAqafCHnvAO+/AHXdAs2aRvmWUawoTgPuAxys7wcwaAqOA6RHmyBg5OTkqBiL13datMHcuHHww9OgBf/lLuAdhl11S8vaRjRTc/XWguotmLweeB1ZFlUNEJGMUFsLRR4e7kP/5z3DsxBNTVhAgxjUFM2sHnAI8FFcGEZG0sGVLWEDu3j1MGT38MOy3XyxR4rwk9W7gGncvMbMqTzSzocBQgI4dO6YgmohIihQXw+GHw5w54RLTe++FvfeOLY75tiZKUby4WTYwxd0PrOC5D4Ft1aANsB4Y6u4vVfWaubm5rtbZIpLxNm/+T2+iu++GTp3glFMiezszm+vuudWdF9v0kbt3dvdsd88GngMura4giIjUCX/7G+y/P0ydGh5fcUWkBaEmorwk9WngLWA/MysysyFmdrGZXRzVe4qIpLUvvoBzz4Xjjw+LxxHcfLazIltTcPezanDuuVHlEBFJCy+8ABdfDGvXwvXXh/bWTZrEnWo76n0kIpIKX30F2dkwY0a4yihNqc2FiEgUSkrgvvvgj38Mj887L+x5kMYFAVQURERq36JFkJcHl1/+n8VkM2jYMN5cSVBREBGpLZs2wU03Qc+e4e7kxx+HiRPjTlUjKgoiIrVl9my4+eaw38GSJXDOOWGEkEG00CwisjO+/RZeew1+9jM46qiwTWaarxtURSMFEZEd9fLLcMABcPrp8O9/h2MZXBBARUFEpOZWrYKzzgodTFu0CNtkxtivqDZp+khEpCbWr4eDDoIvv4QRI+Caa/7Tw6gOUFEQEUnGqlVhB7SsLLjttrAJTrducaeqdZo+EhGpSnEx3HUXdO4c1hAABg+ukwUBNFIQEancu+/CBRdAQUFYP8jwReRkaKQgIlKR226DH/8YVqyAZ5+FyZOhQ4e4U0VORUFEpCJ77hluPluyJNyMlmE3oe0oFQUREQhdTIcOhQcfDI/POy80s2vdOt5cKaaiICLywgth4XjcuHCVUT2mhWYRqb9WroRhw+DFF6FHD5gyBXr1ijtVrDRSEJH6a8kSeOUVGDUK3nmn3hcE0EhBROqbwkLIz4chQ+CYY+Cjj8JNaQJEOFIws/FmtsrM3q/k+V+Y2cLSf2+a2UFRZRERYcsWuPXWcK/BtdfCN9+E4yoI3xPl9NEEoF8Vz38IHOXu3YFbgLERZhGR+qygAHJz4Xe/Cy2uFy4MjexkO5FNH7n762aWXcXzbyY8nA20jyqLiNRja9bAkUdCq1ZhQfnkk+NOlNbSZU1hCDAt7hAiUocsXBimitq0gUmT4IgjoGXLuFOlvdivPjKzvoSicE0V5ww1swIzK1i9enXqwolI5vnii9Cw7qCDYPr0cOy//1sFIUmxFgUz6w48Cgxw9y8qO8/dx7p7rrvntm3bNnUBRSRzuMPTT0PXrvDUU3D99WF7TKmR2KaPzKwj8AJwjrsXxpVDROqIwYPhiSfCPgczZtSLjqZRiKwomNnTQB+gjZkVATcCuwC4+0PAcGB34AELjaaK3T03qjwiUgeVlIRGdQ0awE9/Gm4+u/xyaNgw7mQZy9w97gw1kpub6wUFBXHHEJG4LVoEF14IZ58Nl14ad5q0Z2Zzk/nDO/aFZhGRGtm0CW66CXr2DHcn77573InqlHS5JFVEpHrvvBNaWi9eDL/4BYweDbr4pFapKIhI5vjuO1i3DqZOhf79405TJ6koiEh6mzIldDP97W/h6KPDlFHjxnGnqrO0piAi6enzz2HgwNCr6MknYfPmcFwFIVIqCiKSXtxhwoRwE9qLL8Itt4S1BBWDlND0kYikl48+Cnsl9+4NjzwC++8fd6J6RSMFEYlfcTFMnhw+7twZ3nwTZs1SQYiBioKIxGvBAjj0UBgwAGbPDsdyc8NdypJy+qqLSDw2bIDrrgsF4JNPYOLEMGUksdKagoiknnu4vHT2bDj/fLjjDmjdOu5UgoqCiKTS119D8+Zhaujqq8OWmMccE3cqSaDpIxFJjRdeCJeZPvhgeHzKKSoIaUhFQUSitXIlnHoqnHYa7LlnWFSWtKWiICLRmTgxjA6mTYNRo8JNaD/+cdyppApaUxCR6LRtG4rAww9Dly5xp5EkqCiISO3ZsiVcSbRhQ2hP0bcv9OkTdkeTjKDpIxGpHXPmhHsOrr8eli0Ll52CCkKGUVEQkZ2zbh1cdVVYQF6zBl56CZ5+WsUgQ6koiMjO+eQTeOABuOiisCPagAFxJ5KdEFlRMLPxZrbKzN6v5HkzszFmtszMFppZr6iyiEgt++KLUAggNK1btiw8btky3lyy06IcKUwA+lXxfH+gS+m/ocCDEWYRkdrgHqaGunaFX/867IIG0K5dvLmk1kRWFNz9deDLKk4ZADzuwWxgNzPbO6o8IrKTVqyAE0+EQYNCe+t58yAnJ+5UUsvivCS1HfBJwuOi0mP/Ln+imQ0ljCbo2LFjSsKJSILiYjjySFi9GkaPhssvh4YN404lEYizKFR0aYJXdKK7jwXGAuTm5lZ4johE4F//gn32gUaN4NFHYd99ITs77lQSoTivPioCOiQ8bg+sjCmLiCTatAluvBEOOCDcjQxw7LEqCPVAnEVhMvDL0quQDgW+dvftpo5EJMXefBN69oQRI+DMM+H00+NOJCkU2fSRmT0N9AHamFkRcCOwC4C7PwRMBU4AlgHrgfOiyiIiSbr11nBHcocOMHUq9O8fdyJJsciKgrufVc3zDlwW1fuLSA1s3Ro2vjn00LCIPHIk7Lpr3KkkBmqIJ1KfrVoV7jdo3z40suvbN/yTekttLkTqI3d47LFwE9oLL2h/ZCmjkYJIffPRR3DhhTBjBhxxBDzySGhVIYJGCiL1z8aN8O67oVfRrFkqCPI9GimI1AcLFoRpohEjQhH4+GNo2jTuVJKGNFIQqcs2bIDrrgub34wdC599Fo6rIEglVBRE6qrXXoPu3eG222Dw4LDXwV57xZ1K0pymj0TqonXrwp3Iu+0WFpSPOSbuRJIhNFIQqSvc4a9/hZISaNYMXnkFFi5UQZAaUVEQqQs+/RROPRV++lN46qlwLDcXsrLizSUZR0VBJJNt3Rq6mHbrFkYGo0bBWVV2mBGpktYURDLZBRfAH/8YWlOMHRv2OxDZCSoKIplmy5awbtCkCZx7brgr+bzzwCrat0qkZjR9JJJJ5swJawU33BAeH3kknH++CoLUGhUFkUywbh1ceWVobb1mDeTlxZ1I6ihNH4mku7fegkGDQiO7iy8ON6O1bBl3KqmjVBRE0l2rVqEIvP66RggSOU0fiaQbd3j6abjkkvB4//1h/nwVBEkJFQWRdLJiBZx4Ypgumj8fvv02HNdCsqRIpEXBzPqZ2VIzW2Zm11bwfEczm2lm881soZmdEGUekbRVUgL33htuQps1C+6+G954A5o3jzuZ1DNJrymY2WFAduLnuPvjVZzfELgfOA4oAuaY2WR3X5xw2g3ARHd/0My6AVNL30Okflm7Fm6+Odxz8NBDkJ0ddyKpp5IqCmb2BLAPsAAoKT3sQKVFATgEWObuy0tf4xlgAJBYFBxoUfpxS2Bl0slFMt2mTTB+PAwdCm3aQEEBdOqkqSKJVbIjhVygm7t7DV67HfBJwuMioHe5c24C/mpmlwPNgGNr8PoimeuNN8I+yUuWQOfO0K+fRgeSFpJdU3gfqOnuHBX9uVO+qJwFTHD39sAJwBNmtl0mMxtqZgVmVrB69eoaxhBJI998A5ddFqaJ1q+HadNCQRBJE1WOFMzsL4Rf5M2BxWb2DrBp2/PuflIVn14EdEh43J7tp4eGAP1KX+stM2sCtAFWJZ7k7mOBsQC5ubk1Ga2IpJdTToGZM+GKK+CWW2DXXeNOJPI91U0f3bkTrz0H6GJmnYFPgYHAoHLnrACOASaYWVegCaChgNQtn38efvk3awYjR4Y1g97lZ1JF0kOV00fuPsvdZwEnbPs48Vg1n1sMDAOmA0sIVxktMrMRZrZthHEVcKGZvQs8DZxbw3ULkfTlHtpad+0KN94Yjh16qAqCpDVL5newmc1z917lji109+6RJatEbm6uFxQUpPptRWrmgw/goovg1VfD+sEjj4Q7k0ViYmZz3T23uvOqW1O4BLgU+KGZLUx4qjnw5s5FFKmjJk2CwYOhUSN48MFwyWkDNQ+QzFDdmsJTwDTgViDxjuRv3f3LyFKJZCL3sF7QvXtoVTF6NLRrF3cqkRqpbk3ha3f/yN3PcvePgQ2Eq5F2NbOOKUkoku42bIBrr/3P3sj77QcTJ6ogSEZKakxrZj8zs38BHwKzgI8IIwiR+m3mzDAyGDUqXF20ZUvciUR2SrITnf8HHAoUuntnwmWkb0SWSiTdffVVuCP56KPDtNGMGTBuHOyyS9zJRHZKskVhi7t/ATQwswbuPhPoEWEukfRWXAxTpsDVV8PChXDMMXEnEqkVyfY++srMdgXygSfNbBVQHF0skTT06aehvfXIkaGBXWGhWltLnZPsSGEAsB64AngF+AD4WVShRNLK1q3w8MNhr4N77oEFC8JxFQSpg5IaKbj7OjPrBHRx98fMLAtoGG00kTSwdGlYO8jPD+sHDz8M++4bdyqRyCR79dGFwHPAw6WH2gEvRRVKJC24w5lnwnvvhX0PZsxQQZA6L9k1hcsIm+a8DeDu/zKzPSJLJRKngoLQkmLXXeGJJ2CPPWCvmnaOF8lMya4pbHL3zdsemFkjtt8bQSSzrVsHV14ZGtbdfns41r27CoLUK8mOFGaZ2e+ApmZ2HKEf0l+iiyWSYtOnw8UXw0cfwaWXwm9+E3cikVgkO1K4lrDPwXvARcBU4IaoQomk1B13hN3PmjQJC8r33w8tWlT/eSJ1ULJXH201s5eAl9xdm+BI5nMPPYuysmDAAPjuO7juulAYROqxKkcKFtxkZmuAfwJLzWy1mQ1PTTyRCHz8MZxwApxzTnickwM336yCIEL100dXAIcDB7v77u7eGugNHG5m/xN5OpHaVFICY8bAAQeEaaIjjwwjBhEpU9300S+B49x9zbYD7r7czM4G/gqMjjKcSK1ZvhwGDYK33w7rBw8+CNnZcacSSTvVFYVdEgvCNu6+2szUDlIyR4sW8M038OSTYd8Ds7gTiaSl6qaPNu/gcyLxe+ONsG5QUhIa2L3/fhgtqCCIVKq6onCQmX1Twb9vgR9V9+Jm1s/MlprZMjO7tpJzzjSzxWa2yMye2pH/CZHv+eYbuOwyOOKIsHawYkU4rn2SRapV5fSRu+9w0zszawjcDxwHFAFzzGyyuy9OOKcLcB1wuLuvVesM2WmTJ4ebz1auhCuugFtuCe0qRCQpyd7RvCMOAZa5+3IAM3uG0IJ7ccI5FwL3u/taAHdfFWEeqeu2bAmb3rRqBc8/H9pViEiNRDmebgd8kvC4qPRYohwgx8zeMLPZZtavohcys6FmVmBmBatX6945SeAeFo+/+y5shTltGsydq4IgsoOiLAoVreaVvyi8EdAF6AOcBTxqZrtt90nuY909191z27ZtW+tBJUN98AEcdxycfXbYHxmgc2do3DjeXCIZLMqiUAR0SHjcHlhZwTl/dvct7v4hsJRQJEQqV1wMd94JP/oRvPNOuOfg8svjTiVSJ0RZFOYAXcyss5k1BgYCk8ud8xLQF8DM2hCmk5ZHmEnqgl/9Cn77Wzj+eFiyJHQ31ZVFIrUisoVmdy82s2HAdMLWnePdfZGZjQAK3H1y6XPHm9lioAT4rbt/EVUmyWAbNsD69bD77uGqoqOPhtNO0z0HIrXMPMN6v+Tm5npBQUHcMSSVZs4M+yR37w4vvBB3GpGMZGZz3T23uvM05pb0tXYtXHBBGBUADBsWbx6ReiDK+xREdtzs2XDyybBmDVxzDdx4IzRtGncqkTpPRUHSi3tYJ9h3X+jRA269FXr2jDuVSL2h6SNJD1u3wkMPhfsOtjWwe+UVFQSRFFNRkPgtXQp9+sAll4SRwldfxZ1IpN5SUZD4bNkCI0eGq4refx/Gj4cZM8JlpyISC60pSHy2boU//SksKN9zD+y1V9yJROo9jRQktb77LlxJ9N138IMfwFtvwbPPqiCIpAkVBUmdV16BAw+EESNCN1OA3bbrfygiMVJRkOitWRO2xezfP9xrkJ8PZ5wRdyoRqYCKgkTvoovCFNH//i/Mnx+2yRSRtKSFZonGxx+HfQ323htuvx1uvjlMHYlIWtNIQWpXSQmMGQMHHABXXhmO7bOPCoJIhtBIQWrP+++HBnZvvx3WD267Le5EIlJDGilI7XjpJejVK2yR+eST8PLL0KlT3KlEpIZUFGTnbNoU/puXF/Y8WLIEBg3S5jciGUpFQXbMN9/ApZeGYlBcHFpT3H9/aGQnIhlLRUFqbvJk6NYtdDU9/PBQFESkTlBRkOStXQtnngkDBkDr1qFFxejR0KRJ3MlEpJZEWhTMrJ+ZLTWzZWZ2bRXnnW5mbmbV7h8qMcrKgn/9K3Q2nTsXeveOO5GI1LLIioKZNQTuB/oD3YCzzKxbBec1B34FvB1VFtkJH3wQWlRsa2A3Zw787newyy5xJxORCEQ5UjgEWObuy919M/AMMKCC824Bbgc2RphFaqq4GO64A370o7CGsGBBON5It7aI1GVRFoV2wCcJj4tKj5Uxs55AB3efEmEOqan58+GQQ+Dqq+H442HxYvUrEqknovyzr6IL1b3sSbMGwGjg3GpfyGwoMBSgY8eOtRRPKnXNNbByJUyaBKedpnsOROqRKItCEdAh4XF7YGXC4+bAgcBrFn7p7AVMNrOT3L0g8YXcfSwwFiA3N9eR2jdzJnTpAu3bh20xmzWDVq3iTiUiKRbl9NEcoIuZdTazxsBAYPK2J939a3dv4+7Z7p4NzAa2KwgSsbVrQ7+io4+G3/8+HGvfXgVBpJ6KrCi4ezEwDJgOLAEmuvsiMxthZidF9b6SJHd47jno2hUmTAhTRn/4Q9ypRCRmkV5K4u5Tganljg2v5Nw+UWaRcu67D371q9DEbto06Nkz7kQikgZ0fWF9snUrfPEFtG0bmtaVlMCwYbrMVETKqM1FffHPf8JRR4V9DrY1sLviChUEEfkeFYW6bvNm+L//g4MOgkWLwsigYcO4U4lImtKfiXXZxx/Dz34G770HP/853HMP7Lln3KlEJI1ppFCX7bUX7L13aFPxzDMqCCJSLRWFuuaVV+DII+Hbb0MDu+nTw2hBRCQJKgp1xZo1cPbZYSF59erQpkJEpIZUFDKdOzz5ZLgJbeJEGD48dDTdb7+4k4lIBtJCc13wxz/CPvvAo4/CgQfGnUZEMpiKQiYqKYEHHoBTTgl9iiZOhJYtdampiOw0TR/FpLCwkHHjxlFYWFizT3zvPTjssNCi4rHHwrHWrVUQRKRWaKQQg8LCQnr16oW7Y2bMmzePnJycCs/Lz88nLy+PnI4dw97It90Gu+0GTz0FAwfGkF5E6jIVhRjk5+fj7qxfv56srCzy8/PLjufl5ZGTk7Nd4Vjxi1/QeuzYsF/yXXdBmzYx/1+ISF2kohCDvLw8zIysrCzMjPbt25cVAHfn+uuvB2DXrVvJ2rCBz7OyeLlbN86ZPj1sjykiEhEVhRjk5OQwb968spFB4sgBYMSIEZwEzN28mX83aEAfoHf//lDBFJOISG1SUYhJTk7O99YRzIzGjRvTavNmxmzezJnA53vuyadDhjBv8OAK1xxERGqbrj5KocquONo2cnjgggtYDJwE3LTLLnz96qsMGDlSBUFEUkYjhRSp8oqjLVvCyOGuu/j688+Z0rMng844Q8VARFJORSFFKrriKOeHPwxXEo0bBwUF0Lw5LZ97jtPjDisi9Vak00dm1s/MlprZMjO7toLnrzSzxWa20MxeNbNOUeaJU/krjo5t3RoOOQSuuQa6dYONG+OOKCIS3UjBzBoC9wPHAUXAHDOb7O6LE06bD+S6+3ozuwS4Hfh5VJnikHgD2rx58/jHzJmcXFBA6zPOCHslP/88nHpq3DFFRIBop48OAZa5+3IAM3sGGACUFQV3n5lw/mzg7AjzpFxF6wjnDx0KkybB+efDqFHQqlXcMUVEykQ5fdQO+CThcVHpscoMAaZFmCfltq0j/GD9eu7atIl5f/4zmMHUqTB2rAqCiKSdKEcKVsExr/BEs7OBXOCoSp4fCgwF6NixY23li1zeEUdwakkJdwK7FxezetsTjRvHmEpEpHJRFoUioEPC4/bAdtuBmdmxwPXAUe6+qaIXcvexwFiA3NzcCgtLuti2htC3Sxdy/vAHnti0idWdOvHpPffQacCAuOOJiFQpyqIwB+hiZp2BT4GBwKDEE8ysJ/Aw0M/dV0WYpVZ9r3tpwr0EiWsId2/eTHajRjS4807a/vrX0EhX/4pI+ovsN5W7F5vZMGA60BAY7+6LzGwEUODuk4E7gF2BSWYGsMLdT4oqU20ov3j8/PPPU1RURF5eHgsnTmT/khLmbtzITU2b0nL4cM686qq4I4uIJC3SP1/dfSowtdyx4QkfHxvl+0ch8Sa0pk2bctJJJ9G0YUOuLC7mBmD3khJOzMriazN66FJTEckwmtOoocSb0EpKSugN3L9hAz8Clh98MB3vvpsxS5ZsN7UkIpIJVBRqKLHtdY/PPqPnDTew0owzGjdm5J/+RE5ODvscdljcMUVEdoiKwg7IadOGnCFDYNMm1nz+OX/PyWHk8cdrZCAiGU9FoSbWrIErroBZs2DRImjRgjZjxvDLuHOJiNQS7aeQDHd48kno2hUmToQhQ+AHP4g7lYhIrdNIoQqFhYXM/tvfOO3ZZ2mWnw+HHgqPPgoHHBB3NBGRSKgoVKCwsJBJkyYxcuRIDGi1eTO9b7iBPW66CRo2jDueiEhkVBTKKSws5JwePbh140Z2d6cIGNi0KWOysxmigiAidZyKQqKNG1l31VX8Y8MG1gL7AKsaN8YaNCAvLy/udCIikdNC8zb5+dCjBz2nTGFio0b8uGlT3mnalOHDh39/P2URkTpMI4VtJkyATZtg+nQOzs7mpgoa3omI1HX1uyj8+c/QsSP07AmjR4dF5GbNyAEVAxGpl+rn9NFnn8EZZ8DJJ8Ndd4VjLVpAs2bx5hIRiVn9KgruMG5cuAntL3+B3/8exo+PO5WISNqoX9NHjz0GF1wARx4JjzwCmiISEfmeul8Uiovhww+hSxcYNCjsjzxwIDSoX4MkEZFk1O3fjPPmwSGHQN++sH59KAiDBqkgiIhUom7+dly/Hq6+OhSEf/8bxoyBrKy4U4mIpL26N3302Wdw+OGwfDlceCGMGgWtWsWdSkQkI9SdolBSEu4z2HNPOPbYsG7Qt2/cqUREMkqk00dm1s/MlprZMjO7toLnf2Bmz5Y+/7aZZdf4TdzDHgf77w8rVoAZPPywCoKIyA6IrCiYWUPgfqA/0A04y8y6lTttCLDW3fcFRgOjavQmRUUwYAD8/OfQsmVYSxARkR0W5UjhEGCZuy93983AM8CAcucMAB4r/fg54Bgzs6Re/aGHoFs3mDED7rwTZs8OowUREdlhUa4ptAM+SXhcBPSu7Bx3Lzazr4HdgTXVvnpBAfTuHaaKfvjD2kksIlLPRVkUKvqL33fgHMxsKDAUoGPHjuHgffeFfZKTHFiIiEj1opw+KgI6JDxuD6ys7BwzawS0BL4s/0LuPtbdc909t23btuFgkyYqCCIitSzKojAH6GJmnc2sMTAQmFzunMnA4NKPTwf+7u7bjRRERCQ1Ips+Kl0jGAZMBxoC4919kZmNAArcfTIwDnjCzJYRRggDo8ojIiLVs0z7w9zMVgMfA21IZkE6vWRaZuWNXqZlVt7oRZW5k7u3re6kjCsK25hZgbvnxp2jJjIts/JGL9MyK2/04s5cNxviiYjIDlFREBGRMplcFMbGHWAHZFpm5Y1epmVW3ujFmjlj1xRERKT2ZfJIQUREalnaF4WUtN+uRUnkvdLMFpvZQjN71cw6xZGzXKYqMyecd7qZuZnFejVHMnnN7MzSr/MiM3sq1RkryFPd90VHM5tpZvNLvzdOiCNnaZbxZrbKzN6v5HkzszGl/y8LzaxXqjNWkKm6zL8ozbrQzN40s4NSnbFcnirzJpx3sJmVmNnpqcqGu6ftP8JNbx8APwQaA+8C3cqdcynwUOnHA4Fn0zxvXyCr9ONL4sybbObS85oDrwOzgdx0zgt0AeYDrUof75HuX2PCPPIlpR93Az5xZKblAAAEPUlEQVSKMe+RQC/g/UqePwGYRuhddijwdpxf3yQzH5bw/dA/7szV5U34vvk7MBU4PVXZ0n2kEG377dpXbV53n+nu2zZ+mE3oCRWnZL7GALcAtwMbUxmuAsnkvRC4393XArj7qhRnLC+ZzA60KP24Jdv3CUsZd3+dCnqQJRgAPO7BbGA3M9s7NekqVl1md39z2/cDafBzl8TXGOBy4Hkgpd+/6V4UKmq/3a6yc9y9GNjWfjsOyeRNNITwF1ecqs1sZj2BDu4+JZXBKpHM1zgHyDGzN8xstpn1S1m6iiWT+SbgbDMrIvxleHlqou2Qmn6fp5t0+Lmrkpm1A04BHkr1e6f7Hs211n47RZLOYmZnA7nAUZEmql6Vmc2sAWFXvHNTFagayXyNGxGmkPoQ/iLMN7MD3f2riLNVJpnMZwET3P0PZvYTQk+wA919a/TxaiydfuZqxMz6EorCEXFnqcbdwDXuXpLqiY90Lwo1ab9dVFX77RRJJi9mdixwPXCUu29KUbbKVJe5OXAg8FrpN+dewGQzO8ndC1KW8j+S/Z6Y7e5bgA/NbCmhSMxJTcTtJJN5CNAPwN3fMrMmhB44cU99VSSp7/N0Y2bdgUeB/u7+Rdx5qpELPFP6M9cGOMHMit39pcjfOc7FliQWYxoBy4HO/GeB7oBy51zG9xeaJ6Z53p6ERccucX99k81c7vzXiHehOZmvcT/gsdKP2xCmOnZP88zTgHNLP+5K+CVrMWbOpvJF2//m+wvN78SVswaZOwLLgMPizplM3nLnTSCFC81pPVLwDGu/nWTeO4BdgUmlfwWscPeT0jxz2kgy73TgeDNbDJQAv/UY/zJMMvNVwCNm9j+EqZhzvfQ3QqqZ2dOEqbc2pWscNwK7ALj7Q4Q1jxMIv2TXA+fFkTNREpmHE9YaHyj9uSv2GJvOJZE3NrqjWUREyqT71UciIpJCKgoiIlJGRUFERMqoKIiISBkVBRERKaOiIFKN0i6VC8zsXTObZ2aHlR7/LzN7rpLPyTazQQmPzzWz+1KVWWRHqSiIVG+Du/dw94OA64BbAdx9pbtv19K49M76bGBQ+edE0l1a37wmkoZaAGshjAaAKe5+oJmdS7jTtwnQDMgCuprZAkIX37XAf5nZK8A+wIvufnXK04tUQ0VBpHpNS3+5NwH2Bo6u5LyfAN3d/Usz6wP8xt1PhDB9BPQgtDnZBCw1s3vd/ZNKXkskFpo+Eqnetumj/Ql9lR6vZM+Ov7l7Vc0YX3X3r919I7AYiH3XPZHyVBREasDd3yI02WtbwdPrqvn0xI64JWikLmlIRUGkBsxsf0JTu+oa7H1LaDsuklH0l4pI9batKUBoFz3Yq9/8ZCFQbGbvElofr63qZJF0oS6pIiJSRtNHIiJSRkVBRETKqCiIiEgZFQURESmjoiAiImVUFEREpIyKgoiIlFFREBGRMv8Pxny5SjNDU/MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hom_dim = 1\n",
    "pairs = st.persistence_pairs()\n",
    "simps = [s for s in filter(lambda x: len(x[0]) == hom_dim + 1,pairs)]\n",
    "pd = np.array([[st.filtration(x),st.filtration(y)] for x,y in simps])\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(pd[:,0],pd[:,1], s = 8 ,color = 'black')\n",
    "plt.plot([0,1.5],[0,1.5],ls = '--',color = 'red')\n",
    "plt.xlim([-0.01,1.5])\n",
    "plt.ylim([-0.03,1.5])\n",
    "plt.xlabel('Birth')\n",
    "plt.ylabel('Death')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bayes_tda module works with persistence diagrams given in birth-persistence (where persistence = death - birth) instead of birth-death coordnates. We can convert our diagram to birth-persistence coordinates in one line of code then visualize the transformed diagram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFl9JREFUeJzt3X+U5XV93/Hni4UNroJadzwYlnVpy0SJQaETopiNUE2y0ByoLSSgpCFnj3uMv3qqtVHTqtGjpKZWS4JFKhzUnkBAG7O1WNJjNExVDMOqVNawZ4s/mFLKikRbF8Rd3/3j3vlmdpiZ+511v/c7yz4f5+zZ+733c++8ds7sfc33+7nfzzdVhSRJAEf1HUCStHpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWoc3XeAlVq/fn1t2rSp7xiSdFi54447vl1VE6PGHXalsGnTJmZmZvqOIUmHlSTfbDPOw0eSpIalIElqWAqSpIalIElqWAqSpIalIElqWAqSpIalIElqWAqSpEZnpZDk2iQPJPnqiHE/m2R/kgu7yiJJaqfLPYXrgC3LDUiyBvjXwC0d5pAktdRZKVTVrcB3Rgx7LfBx4IGuckiS2uttTiHJicBLgav6yiBJOlCfE83vB367qvaPGphkW5KZJDN79uwZQzRJOjL1uXT2FHBDEoD1wHlJ9lXVJxYOrKqrgasBpqamaqwpJekI0lspVNXJc7eTXAd8crFCkCSNT2elkOR64GxgfZJZ4G3AMQBV5TyCJK1CnZVCVV2ygrGXdZVDktSeZzRLkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqWgiSpYSlIkhqdlUKSa5M8kOSrSzz+8iR3Dv98Pslzu8oiSWqnyz2F64Atyzz+deBFVXUa8E7g6g6zSJJaOLqrF66qW5NsWubxz8/bvA3Y0FUWSVI7q2VOYSvwqb5D9G3Xrl1cc8017Nq1q+8oko5Qne0ptJXkHAal8PPLjNkGbAPYuHHjmJKN165duzjjjDOoKpKwY8cOJicn+44l6QjT655CktOADwEXVNWDS42rqquraqqqpiYmJsYXcIymp6epKvbu3UtVMT093XckSUeg3kohyUbgPwG/XlVH/PGSzZs3k4R169aRhM2bN/cdSdIRqLPDR0muB84G1ieZBd4GHANQVVcBbwWeBnwgCcC+qprqKs9qNzk5yY4dO5ienmbz5s0eOpLUi1RV3xlWZGpqqmZmZvqOIUmHlSR3tPnFe7V8+kiStApYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkhqUgSWpYCpKkRmelkOTaJA8k+eoSjyfJFUl2J7kzyRldZZEktdPlnsJ1wJZlHj8XOGX4Zxvw7zvMIklqobNSqKpbge8sM+QC4CM1cBvwlCTP6CqPJGm0PucUTgTunbc9O7zvMZJsSzKTZGbPnj1jCSdJR6I+SyGL3FeLDayqq6tqqqqmJiYmOo4lSUeuPkthFjhp3vYG4L6eskiS6LcUtgP/ZPgppOcD362q/91jHkk64h3d1QsnuR44G1ifZBZ4G3AMQFVdBdwMnAfsBvYCv9lVFklSO52VQlVdMuLxAl7d1deXJK2cZzRLkhqWgiSpYSlIkhqWgiSpYSlIkhorKoUkT+wqiCSpf61KIclZSXYCXxtuPzfJBzpNJkkau7Z7Cu8Dfhl4EKCqvgL8QlehJEn9aH34qKruXXDX/kOcRZLUs7ZnNN+b5CygkqwFXsfwUJIk6fGj7Z7CKxksSXEig9VNn4dLVEjS406rPYWq+jbw8o6zSJJ61vbTRx9O8pR5209Ncm13sSRJfWh7+Oi0qvrruY2qegg4vZtIkqS+tC2Fo5I8dW4jyd+iw2W3JUn9aPvG/l7g80k+Nty+CHhXN5EkSX1pO9H8kSR3AOcAAf5RVe3sNJkkaexWcgjor4CH5p6TZGNVfauTVJKkXrQqhSSvZXCN5f/D4EzmAAWc1l00SdK4tZ1o/qfAT1XVT1fVaVX1M1U1shCSbElyd5LdSd60yOMbk3wmyZeS3JnkvJX+AyRJh07bUrgX+O5KXjjJGuBK4FzgVOCSJKcuGPYvgRur6nTgYsCVVyWpR23nFO4BPpvkvwA/mLuzqv7tMs85E9hdVfcAJLkBuACYP0FdwPHD208G7muZR5LUgbal8K3hn7XDP22cyGAPY84s8HMLxrwd+LPhnMUTgZe0fG1JUgfafiT1d2Fw5bWq+n7L185iL7Vg+xLguqp6b5IXAB9N8pyq+tEBL5RsA7YBbNy4seWXlyStVNu1j15wEFdemwVOmre9gcceHtoK3AhQVV8AjgXWL3yhqrq6qqaqampiYqJNZEnSQWg70fx+Vn7ltduBU5KcPLwGw8XA9gVjvgW8GCDJsxmUwp6WmSRJh1hnV16rqn3Aa4BbGOxh3FhVdyV5R5Lzh8PeALwiyVeA64HLqmrhISZJ0ph0euW1qroZuHnBfW+dd3sn8ML2cSVJXfpxrrz2qq5CSZL60XZP4aeq6oArryV5IfC5Qx9JktSXtnsKf9DyPknSYWzZPYXhuQNnARNJXj/voeOBNV0GkySN36jDR2uBJw3HHTfv/u8BF3YVSpLUj2VLoar+AviLJNdV1TcBkhwFPKmqvjeOgJKk8Wk7p3B5kuOTPJHBgnZ3J3ljh7kkST1oWwqnDvcM/iGD8w42Ar/eWSpJUi/alsIxSY5hUAp/WlU/5LGL20mSDnNtS+GDwDcYLG99a5JnMphsliQ9jrRdOvsK4Ip5d30zyTndRJIk9WXUeQqXVtV/XHCOwnzLXXlNknSYGbWn8MTh38ctO0qS9Lgw6jyFDw7//t3xxJEk9antldfeMzxP4Zgkn07y7SSXdh1OkjRebT999EvD8xR+hcHS2ZOAJ69J0uNM6/MUhn+fB1xfVd/pKI8kqUdtr6ewPclfAQ8Dr0oyATzSXSxJUh9G7ikMF8D7z8ALgKnh2cx7gQs6ziZJGrORpVBVPwLeW1UPVdX+4X3fr6r7O08nSRqrtnMKf5bkHyfJSl48yZYkdyfZneRNS4z51SQ7k9yV5I9W8vqSpEOr7ZzC6xmcyLY/ycNAgKqq45d6QpI1wJXALzL4xNLtSbZX1c55Y04B3gy8sKoeSvL0g/x3SJIOgbZrHx3MGc1nArur6h6AJDcwmIfYOW/MK4Arq+qh4dd54CC+jiTpEGl78lqSXJrkXw23T0py5oinnQjcO297dnjffJPAZJLPJbktyZYlvv62JDNJZvbs2dMmsiTpILSdU/gAg08fvWy4/f8YHBpazmLzDwuvwXA0cApwNnAJ8KEkT3nMk6qurqqpqpqamJhoGVmStFJtS+HnqurVDM9NGB7uWTviObPASfO2NwD3LTLmT6vqh1X1deBuBiUhSepB21L44XDiuACGJ6/9aMRzbgdOSXJykrXAxcD2BWM+AZwzfM31DA4n3dMykyTpEGtbClcAfwI8Pcm7gP8OvHu5J1TVPuA1wC3A14Abq+quJO9Icv5w2C3Ag0l2Ap8B3lhVDx7Ev0OSdAikqt2llpM8C3gxg7mCT1fV17oMtpSpqamamZnp40tL0mEryR1VNTVq3Kgrrx0LvBL4u8D/AD443AOQJD0OjTp89GFgikEhnAv8m84TSZJ6M+rktVOr6mcAklwD/GX3kSRJfRm1p/DDuRseNpKkx79RewrPTfK94e0ATxhuj1z7SJJ0+Fm2FKpqzbiCSJL61/Y8BUnSEcBSkCQ1LAVJUsNSkCQ1LAVJUsNSkCQ1LAVJUsNSkCQ1LAVJUsNSkCQ1LAVJUsNSkCQ1LAVJUqPTUkiyJcndSXYnedMy4y5MUklGXj9UktSdzkohyRrgSgaX8TwVuCTJqYuMOw54HfDFrrJIktrpck/hTGB3Vd1TVY8CNwAXLDLuncB7gEc6zCJJaqHLUjgRuHfe9uzwvkaS04GTquqTHeaQJLXUZSlkkfuqeTA5Cngf8IaRL5RsSzKTZGbPnj2HMKIkab4uS2EWOGne9gbgvnnbxwHPAT6b5BvA84Hti002V9XVVTVVVVMTExMdRpakI1uXpXA7cEqSk5OsBS4Gts89WFXfrar1VbWpqjYBtwHnV9VMh5kkScvorBSqah/wGuAW4GvAjVV1V5J3JDm/q68rSTp4R3f54lV1M3DzgvveusTYs7vMIkkazTOaJUkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEkNS0GS1LAUJEmNTkshyZYkdyfZneRNizz++iQ7k9yZ5NNJntllHknS8jorhSRrgCuBc4FTgUuSnLpg2JeAqao6DfgY8J6u8kiSRutyT+FMYHdV3VNVjwI3ABfMH1BVn6mqvcPN24ANHeaRJI3QZSmcCNw7b3t2eN9StgKf6jCPJGmEozt87SxyXy06MLkUmAJetMTj24BtABs3bjxU+SRJC3S5pzALnDRvewNw38JBSV4C/A5wflX9YLEXqqqrq2qqqqYmJiY6CStJ6rYUbgdOSXJykrXAxcD2+QOSnA58kEEhPNBhFklSC52VQlXtA14D3AJ8Dbixqu5K8o4k5w+H/T7wJOCmJF9Osn2Jl5MkjUGXcwpU1c3AzQvue+u82y/p8utLklbGM5olSQ1LQZLUsBQkSQ1LQZLUsBQkSQ1LQZLUsBQkSQ1LQZLUsBQkSQ1LQZLUsBQkSQ1LQZLUsBQkSQ1LoSe7du3immuuYdeuXX1HkaRGp0tna3G7du3ijDPOoKpIwo4dO5icnOw7liS5p9CH6elpqoq9e/dSVUxPT/cdSZIAS6EXmzdvJgnr1q0jCZs3b+47kiQBlsJYzc0jAOzYsYMrrriiOXTkHIOk1cA5hYOwa9cupqen2bx5c+u5gMXmEbZu3brkY84xSOqDpbBCK30DnyuQ+++/v5lHWLduHdPT083z5s8xLHxMksap01JIsgX4d8Aa4ENV9XsLHv8J4CPA3wMeBH6tqr7RZaYf18I38JtuuokTTjhh0b2G+QVSVQCLziPMn2OoKu6///7mMNJK90gk6ceRuTerQ/7CyRpgF/CLwCxwO3BJVe2cN+ZVwGlV9cokFwMvrapfW+51p6amamZmppPMc+Z+u9+wYQOzs7MHvCkv9kafpNlrgL95I5+enuZ1r3tdUyBvectbli2Qm266iXe9610kWfS1LQZJByvJHVU1NWpcl3sKZwK7q+qeYaAbgAuAnfPGXAC8fXj7Y8AfJkl11VQtzL3p79+/n0ceeYQnPOEJHHXUUc2b8uTkJDt27GgOCb373e8+YK/h8ssvbw4tffzjHz/gU0YXXXTRkm/sk5OTnHDCCSRh7969rF27FoBHH33UQ0qSxqbLTx+dCNw7b3t2eN+iY6pqH/Bd4GkdZhpp7vDQI488AsDDDz/8mHMJJicn2bp1KxdddNEBb/rAAecfzM7OPuZTRsuZfxhpzZo1rFmzxo+tShqrLvcUssh9C/cA2owhyTZgG8DGjRt//GTLmHtjPvbYY5s9haXelOfvNcw9fvnllx/wRj63d9HGYq/nnIKkcepyTuEFwNur6peH228GqKrL5425ZTjmC0mOBu4HJpY7fNT3nELb5/pGLmk1WQ1zCrcDpyQ5GfhfwMXAyxaM2Q78BvAF4ELgz/ucT5izkt/uD+VzJalvnZVCVe1L8hrgFgYfSb22qu5K8g5gpqq2A9cAH02yG/gOg+KQJPWks8NHXUmyB/gmsB74ds9xVupwy2ze7h1umc3bva4yP7OqJkYNOuxKYU6SmTbHx1aTwy2zebt3uGU2b/f6zuyCeJKkhqUgSWoczqVwdd8BDsLhltm83TvcMpu3e71mPmznFCRJh97hvKcgSTrEVn0pJNmS5O4ku5O8aZHHfyLJHw8f/2KSTeNPeUCeUXlfn2RnkjuTfDrJM/vIuSDTspnnjbswSSXp9dMcbfIm+dXh9/muJH807oyL5Bn1c7ExyWeSfGn4s3FeHzmHWa5N8kCSry7xeJJcMfy33JnkjHFnXCTTqMwvH2a9M8nnkzx33BkX5Fk277xxP5tkf5ILx5WtWQJ6Nf5hcNLb/wT+NrAW+Apw6oIxrwKuGt6+GPjjVZ73HGDd8PZv9Zm3bebhuOOAW4HbgKnVnBc4BfgS8NTh9tNX+/eYwXHk3xrePhX4Ro95fwE4A/jqEo+fB3yKwdplzwe+2Of3t2Xms+b9PJzbd+ZReef93Pw5cDNw4biyrfY9hWb57ap6FJhbfnu+C4APD29/DHhx5pYsHb+ReavqM1W1d7h5G7BhzBkXavM9Bngn8B7gkXGGW0SbvK8ArqyqhwCq6oExZ1yoTeYCjh/efjJw3xjzHRik6lYGKwws5QLgIzVwG/CUJM8YT7rFjcpcVZ+f+3lgFfy/a/E9Bngt8HFgrD+/q70UDrflt9vknW8rg9+4+jQyc5LTgZOq6pPjDLaENt/jSWAyyeeS3Da8AmCf2mR+O3BpklkGvxm+djzRDspKf85Xm9Xw/25ZSU4EXgpcNe6vvdqv0XzIlt8ek9ZZklwKTAEv6jTRaMtmTnIU8D7gsnEFGqHN9/hoBoeQzmbwG+F0kudU1V93nG0pbTJfAlxXVe8drjD80WHmH3Ufb8VW0/+5FUlyDoNS+Pm+s4zwfuC3q2r/uA98rPZSmAVOmre9gcfuVs+NmR0uv/1kRu+WdaVNXpK8BPgd4EVV9YMxZVvKqMzHAc8BPjv84TwB2J7k/Krqdg3zxbX9mbitqn4IfD3J3QxK4vbxRHyMNpm3AlsAarCU/LEM1sDp+9DXYlr9nK82SU4DPgScW1UP9p1nhCnghuH/ufXAeUn2VdUnOv/KfU62tJiMORq4BziZv5mg++kFY17NgRPNN67yvKczmHQ8pe/vb9vMC8Z/ln4nmtt8j7cAHx7eXs/gUMfTVnnmTwGXDW8/m8GbbHrMvImlJ23/AQdONP9lXzlXkHkjsBs4q++cbfIuGHcdY5xoXtV7CnWYLb/dMu/vA08Cbhr+FvCtqjp/lWdeNVrmvQX4pSQ7gf3AG6vH3wxbZn4D8B+S/DMGh2Iuq+E7wrgluZ7Bobf1wzmOtwHHAFTVVQzmPM5j8Ca7F/jNPnLO1yLzWxnMNX5g+P9uX/W46FyLvL3xjGZJUmO1f/pIkjRGloIkqWEpSJIaloIkqWEpSJIaloI0wnCVyi8n+UqSHUnOGt7/k0k+tsRzNiV52bzty5L84bgySwfLUpBGe7iqnldVzwXeDFwOUFX3VdVjljQenlm/CXjZwsek1W5Vn7wmrULHAw/BYG8A+GRVPSfJZQzO9D0WeCKwDnh2ki8zWMX3IeAnk/xX4O8Af1JV/2Ls6aURLAVptCcM39yPBZ4B/P0lxr0AOK2qvpPkbOCfV9WvwODwEfA8Bsuc/AC4O8kfVNW9S7yW1AsPH0mjzR0+ehaDdZU+ssQ1O/5bVS23GOOnq+q7VfUIsBPo/ap70kKWgrQCVfUFBovsTSzy8PdHPH3+irj7cU9dq5ClIK1AkmcxWNRu1AJ7/5fBsuPSYcXfVKTR5uYUYLBc9G/U6Iuf3AnsS/IVBksfP7TcYGm1cJVUSVLDw0eSpIalIElqWAqSpIalIElqWAqSpIalIElqWAqSpIalIElq/H8UF+At8czm+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pd_t = np.array([[b,d-b] for b,d in pd])\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(pd_t[:,0],pd_t[:,1], s = 8 ,color = 'black')\n",
    "plt.xlim([-0.01,1.5])\n",
    "plt.ylim([-0.03,1.5])\n",
    "plt.xlabel('Birth')\n",
    "plt.ylabel('Persistence')\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our Bayesian model views a persistence diagram $\\mathcal{D}$ as a sample from a point process $D$, which is a random collection of elements whose cardinality may vary. A useful summary statistic of a point process is its intensity,$\\lambda$, which satisfies $\\mathbb{E}(|D \\cap \\mathcal{A}|) = \\int _{A}\\lambda(x)\\,dx$ for every $A \\subset \\mathbb{W} := \\{(b,p) \\in \\mathbb{R}^2:b >0\\}$. Our Bayesian model gives a recipe to compute posterior intensities using prior intensities, a likelihood, and observed persistence diagrams. In our implementation, the prior intensity $\\lambda _D(x)$ is given by a Gaussian mixture. \n",
    "\\begin{equation}\n",
    "\\lambda\\_D(x) = \\sum\\_{i=1}^{N} w_i\\mathcal{N}^'(x;\\mu_i,\\sigma_i I), \n",
    "\\end{equation}\n",
    "where $N$ is the number of mixture components, $w_i$ are weights, and \n",
    "\\begin{equation}\n",
    "\\mathcal{N}^'(x;\\mu_i,\\sigma_i I) = \\begin{cases}\\mathcal{N}(x;\\mu_i,\\sigma_i I), x \\in \\mathbb{W} \\newline\n",
    "                                                    0, \\mathrm{otherwise}\n",
    "                                    \\end{cases}\n",
    "\\end{equation}\n",
    "where $\\mathcal{N}(x;\\mu_i,\\sigma_i I)$ is the Gaussian density, with mean $\\mu _i$ and covariance matrix $\\sigma _i I$. Thus, we can create prior intensities by specifying the means, covariances, and weights of each mixture component, respectively. Below, we demonstrate the syntax for doing so."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bayes_tda import *\n",
    "\n",
    "mu_pri = np.array([[0.2,1]]) # prior mean\n",
    "w_pri = np.array([1]) # prior weight\n",
    "sig_pri = np.array([0.25]) # prior covariance magnitude\n",
    "\n",
    "prior_unimodal = Prior(weights = w_pri,mus = mu_pri, sigmas = sig_pri)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object prior_unimodal corresponds to the prior intensity\n",
    "\\begin{equation}\n",
    "\\lambda\\_D(x) =  \\mathcal{N}^'\\Big(x;(0.2,1),0.25 I\\Big)\n",
    "\\end{equation}\n",
    "\n",
    "We can evaluate the prior by calling\n",
    "```python\n",
    "prior_unimodal.eval(x)\n",
    "```\n",
    "e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.6366197723675814, 0.17349912804424186, 0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prior_unimodal.eval([0.2,1]), prior_unimodal.eval([0.1,0.2]), prior_unimodal.eval([1,-0.1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the point $(1,-0.1)$ evaluates to 0. Since we work with homologoical features in birth-persistence coordinates, our model assigns probability zero to any point with negative persistence. Another component of our model is an intensity for unanticipated features,which measures the probability that features arise from noise instead of characteristic topology in the data. This intensity, which is another Gaussian mixture supplied by the user, serves as another source of prior information and prevents overfitting the posteriors to noise. In this example since the noise level is relatively low, we anticipate that features whose persistences are close to zero are the result of noise. We leverage this in creating our intensity for unanticipated features.     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_u = np.array([[0.15,0.01],[0.2,0.01],[0.25,0.01]]) # unanticipated features means\n",
    "w_u = np.array([2,2,2]) # unanticipated features weights\n",
    "sig_u = np.array([0.05,0.05,0.05]) # unanticipated features covariance magnitudes\n",
    "\n",
    "prior_unant = Prior(weights = w_u,mus = mu_u, sigmas = sig_u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last ingredient we need to compute a posterior distribution is a likelihood $\\ell$, which in our implementation is given by a restricted Gaussian density, \n",
    "\\begin{equation}\n",
    "\\ell(y|x) = \\mathcal{N}^'(y;x,\\sigma\\_L I)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [],
   "source": [
    "sl = 0.3 # magnitude of covariance matrix for likelihood"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we compute the posterior distribution with one last line of code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "posterior = Posterior(prior = prior_unimodal,clutter = prior_unant, Dy = pd_t, sy = sl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Respectively assuming that our prior and likelihood is a restricted Gaussian mixture and a restricted Gaussian density results in a posterior that is a resticted Gaussian mixture (this is known as conjugacy in Bayesian inference). The exact equations for the posterior means, weights, and covariance matrices may be found in https://epubs.siam.org/doi/pdf/10.1137/19M1268719. Let's visualize our posterior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXu8XGV577+/vXMFQpJNuIarNSCIFzBAUat4o4HTI3pKe0BpxUOlHEXbKm31lA8ilI/W1muLxdRykPYIKj3a1IOiFSlWRBOroATQiFxCqJALCZckZO/9nD/WmrAymcu71qzrzPP9fOazZ2a9a61nZs981zPP+653ycxwHMdxhpexqgNwHMdxisVF7ziOM+S46B3HcYYcF73jOM6Q46J3HMcZclz0juM4Q46L3qkFkp6U9Jyq4ygKSb8m6d6q43CeRdK5kv696jjKwEU/Yki6X9LWWKy/lPS/Je01wPYOl2SSZgwSl5ntZWb3DbKNREyXSvrHwLalfNnN7NtmdlRiv/dLem3R+x0WJL1c0m2SNkvaKOk7kk4ocf/BnxNJp0haW3RMaXDRjyb/1cz2Ao4HTgAuriqQQQ8Qg67v1B9JewNfAf4amAAWAx8AtlcZV5Nw0Y8wZvYw8FXgWABJB0laEWdMayS9rdVW0omSVknaEv8S+Gi86Nb47+Pxr4ST4/b/Q9LdkjZJuknSYYltmaR3SPoZ8LPEc8+N78+XdK2kxyQ9IOliSWPxsnPjbO5jkjYCl/Z7nfG2L5D0szieKxVxNHAVcHIc++Nx+9mS/krSg/FrvUrS3HjZKZLWSnqPpEclPSLprYl9nS5ptaQnJD0s6aLkevH9fwAOBf4l3u+fSPp/kt7ZFvedkt4Q8K8cdo4EMLPrzGzKzLaa2dfN7M5WA0lviz9vT8Tv//Hx8++V9PPE82/MI6D4F9lF8f9os6TPS5ojaU+i79RB8f/2yfh7NZaIZYOkL0iaiLfV+lX8lvgzt17SnyX21e27F46Z+W2EbsD9wGvj+4cAdwGXx4//DfgUMAd4MfAY8Jp42XeB34nv7wX8anz/cMCAGYl9vAFYAxwNzCD6xXBbYrkB3yDKzuYmnntufP9a4J+BefH2fwqcFy87F5gE3hlve26H13gp8I9t+/sKsIBIsI8ByxLb+/e29T8OrIjjmwf8C/DBeNkp8f4vA2YCpwNPAwvj5Y8AvxbfXwgcn1hvbaf/Q/z4t4HvJR6/CNgAzKr6M1P1Ddg7fi8+C5zWeq8Ty38LeJjo16mA5wKHJZYdRJTU/nfgKeDAbv/7HjHs0jb+/30/3vYEcDdwQaf/dfzcHwK3AwcDs4FPA9e1fYf+Dpgb/++3A0f3+u6leg+r/if6rdxb/AF9EngceIBI7HOJpD8FzEu0/SBwTXz/VqKfy4vattf6kCZF/1ViMcePx2IZHhY/NuDVbdux+As6Hn/Ij0ks+33glvj+ucCDfV7jpewu+pcnHn8BeG9ie8kvsGIZ/EriuZOBX8T3TwG2tr3eR3n2wPdgHO/ebTHt8uVnd9HPBjYCS+LHfwV8qurPS11uREnDNcBaogPtCmD/eNlNwB8EbudHwBmd/vd91mv/nNwPnJN4/GHgqk7/6/i5u4mTpvjxgcAOomSl9R06OLH8+8BZ8f2O3700Ny/djCZvMLMFZnaYmb3dzLYSZSYbzeyJRLsHiOqhAOcR/YS+R9JKSb/RY/uHAZ+Q9HhcDtlIJNDFiTYPdVl3ETAr3nenOHqt24v/TNx/migz6sS+wB7ADxLxfy1+vsUGM5vssr3fJMryH5D0b61SVj/MbDvRAeicuEx1NvAPIeuOAmZ2t5mda2YHE5UaDyL65QVRkvLzTutJ+l1JP0r8L48l+ozlQehnCqLvxJcScdxNlFjtH7C9NN+9jrjonRbrgAlJ8xLPHUr0kxgz+5mZnQ3sB/wFcENcj+w0/elDwO/HB5PWba6Z3ZZo023a1PVEmc5hied2xtFn3Sy0b2s9Ucb+/ETs8y3qvO6/MbOVZnYG0fv0ZSJ5h+wXotLEm4HXAE+b2XeDXsGIYWb3EGX3x8ZPPQT8Snu7uF/o74ALgX3MbAHwE6Kko9AQOzz3EHBa23dijkX9ZL031v27F4yL3gHAzB4CbgM+GHcqvZAok/g/AJLOkbSvmU0TlX0gykgeA6aB5Bj4q4D3SXp+vO58Sb8VGMcUkRyvkDQv/rK+GwgaLpmBXwIHS5oV73+aSA4fk7QfgKTFkn6934YkzZL0ZknzzWwHsIXoPeq2313OG4jFPg18BM/mdyLpeXHn98Hx40OIfvHcHjf5DHCRpJfEnezPjT83rUTksXi9t/LswaFIfgnsI2l+4rmriD7Th8Wx7CvpjJCN9fjuBeOid5KcTVQvXAd8CXi/mX0jXrYMuEvSk8AniOqH28zsaeAK4Dvxz9JfNbMvEWUe10vaQpRFnZYijncS1cnvA/4d+Bxw9cCvrjM3E3VI/6ek9fFzf0rUmXx7HP+/Akd1Wb+d3wHuj9e7ADinS7sPAhfH79lFieevBV5AcQe2JvIEcBLwPUlPEQn+J8B7AMzsi0Sfwc/Fbb8MTJjZaqKD5neJ5PsC4DvddiLpLklvHjTY+BfHdcB98f/3IKLvzArg65KeiF/DSYGb7Pjdi2N+UtKv9duA4mK/4zg1QNLvAueb2curjsUZHjyjd5yaIGkP4O3A8qpjcYaLwkQv6WpFJ5T8pE+7EyRNSTqzqFgcp+7EfQCPEZUYPldxOM6QUVjpRtIriMZrX2tmHTtAJI0TnTizDbjazG4oJBjHcZwRprCM3sxuJRo/3Yt3Av9EdMKJ4ziOUwCVTQglaTHwRuDVRKcu92p7PnA+wB5z57zkyCMOfXbZ1GS31bpvL8M6dcTG0//7QteZCpwrbLpPuynrn0v0azM51X/Y89R03yZMpvi3T4VscJdtp2u/2/4mU42WS7ftHcPxec/CE5vuWW9m+/Zv2Z0T586zzdNh/5+fPrP1JjNbNsj+iqDKmf8+DvypmU1Jvb/IZracuIPquOcfZTd/bjmzn1zfc512xjZvyBpnX2xjWCyayOuEvN2Znr9Pqvbb9wqL5am5/bf7xPjCnss37+h/rtGW7XN7Lt/41My+29i0pf8BYePj4ULduHFbcFuAxzc8nap9kk2Pbc68bj82P9rvh/XwcvPnT3qgf6vebJ6e4m8PfG5Q29c88OPivuQDUKXolxKNs4bolOTTJU2a2Zd7raSpyVSSz0PwoSLPYztZDwbJ1xki/dZ72E/4e27d0Ff286Y29ZT9/JlPBsm+FxN77ugr+4V7W1/ZTywYTyX7sli47/zCZD9/v4mRlr1ToejN7IjWfUnXAF/pJ/lQBpF7XlLPa/9ZxN96/XkJf8+t0fZ6CX9Q2e89e2vfrD5E9iGEyn5iYk6qrH7BPnsMlNU7TlEUObzyOqIz0o5SNH/3eYrmBL+gqH2Obd6QSfK2cf3OW91IxpY2vjTvR8ivpJbwuzFvalPP5fNnPtlz+d6zt/aNoR8L9853FNnExJxU7Rfss0fmfS3cd37/Ro6TgcIy+ngSntC25w6yr6xybyLJuEOz/dAMPzS7rzKzH/YSTlF4+Wa0afxl2NJIvii571gftt2Zi/Ltp2m9niKE33TZ50mZJZwia/XO6NJY0Zct+FCZD7KNrAeCtFn+2OYNucgeutftq5a9Z/W741n96NI40WtqMljyWQWfh9Tz2G8W8Ydm+SHZ/aClnDI6aHuRp+y9Y9ZpMkM5qVmWjssd69fvvNWFQWIKfQ9CDpr9Omp7ddIO2kHbi4k9d2Ret854p6yTN0Ml+rSCr6Pcu5E11lDZ9xN+VbLvNxKnn+xDRuFMLBjv2wbSj8CpI/P3m6g6BKcChkL0WQXfVNJKP6/svkjZ92LQYZd5D7kMxYdaOnWhcTX6dtIKflC2/TLfA8Sc/QcbidN6TSH1fNu4Pqh2369uX0TNfpCzZ/MYhVNUrd5pPuOzx5k4MvDAO/CEC8XQWNGXIfi8pR66jyzyDxV+SGftoLLvRVbZlzEKpwjqONTSR9+MHo0s3YRKPkuJZtsv1++8VcUgMYS+5n7v4SBlnH5n0GYljzNnezFKtXpntGic6C1giuG0gq+D3LuRNba8ZN9L+Fll36teX+QonCbW6h0nDxon+n5kEXxTSBtvyAEvpKO2TrIftay+qE5ZH30zWjS2Rt9OWsHnzZOP9K557nVgfl+sZPwh9fwd69cH1e6z1u171ewHOaGqG73q9XWt1TtOlQyF6EMlP6jg+8k867qDHAS2/XJ9abLvxSAdtJ3IYw77rBQxNYKfKetUSeNLNyGSH6RE8+QjG3feiiK5jyz7CX19oaWcbmSd57/sEk5ZtXov3zhpkbRM0r2S1kh6b4flh0n6pqQ7Jd0i6eDEsrdI+ll8e0ua/TZa9KGST0sZci9i/2mE34ussq9iJI7jNAVJ48CVwGnAMcDZko5pa/ZXwLVm9kLgMuCD8boTwPuBk4ATgfdLCq57NlL0IZlp2iy+arl3I6vw+1GF7LtRRVbfj9BO2TT46JuR50RgjZndZ2bPANcDZ7S1OQb4Znz/W4nlvw58w8w2mtkm4BtA8EXIGyd6C7iifRbB1520cVYp+24UNUVCFppWvnEawyJJqxK38xPLFgMPJR6vjZ9Lcgfwm/H9NwLzJO0TuG5XhqIzNklayTeNVswhHbghHbUhnbRpyToKpxtZO2YHnRqhTvPV+wVJqmN81jh7HxTcT7LezJZ2WdZpuFd7xnER8DeSzgVuBR4GJgPX7cpQiT5U8nkJfsu67F+8FB+cjjz5yMZg2UPvYZi9ZN9rJE7WIZfdyDLccpA566sYalm30Tc+HUKprAUOSTw+GFiXbGBm64D/BiBpL+A3zWyzpLXAKW3r3hK648aVbroRIvlByzRb1m3e5TYIeWwrzevp9/70KuOUWcLpxiBnzBaNl2+cQFYCSyQdIWkWcBawItlA0iJJLS+/D7g6vn8TcKqkhXEn7Knxc0EMhehDJZ+FvMRe5H6qln038u6Y7UbTOmWd0cTMJoELiQR9N/AFM7tL0mWSXh83OwW4V9JPgf2BK+J1NwKXEx0sVgKXxc8F0fjSTVGSL1rsIftOU94Jrd2HnmCVhpBr0LaTpVZfBE06U9br9M3HzG4Ebmx77pLE/RuAG7qsezXPZvipKCyjl3S1pEcl/aTL8jfHJwXcKek2SS9Ku48iJF9G9h5KllgG7X+oQ1bfjWEp3/gwS6dsiizdXEPvcZ6/AF4ZnxhwObA8zcbzlnzRgt/4zA7u2vwEG59JX0rIW/aDlHC6kfWs2U7UrXwz7PgZssNPYaI3s1uBrsYxs9vigf8AtxP1IgeRp+TLyOC/+sijnP7tlVzwg59w+rdX8rVHHk29jbRxDir7buSZ1delU7bfmHqv0ztNpy6dsecBX+22UNL5rRMQNjzVf5raNJIvmo3P7OADq9ewfXqaJ6em2D49zaWr12TK7CHfmHvJPksJp+qs3nGczlTeGSvpVUSif3m3Nma2nLi08+ID9+mZflUh+Y0/7b6te7Y/zQyD7YnnZko8snUbE7OyncyzZd3moI7a0LH2ZdBtXH2enbKDTF+cB2muJ5t1PL13yDpZqFT0kl4IfAY4zcwGTgfLknwvsbdzwIxZTLadwLZjepq5656BAc6Zykv2vUbhZDmRKssIHMepM+Mzx2uTMGWlstKNpEOB/wv8jpn9tKz9DiL5jT/dnEryAAvGZ3DRxMHMlthDY8yWuGjiYBaMz9i5vbTbbBH6Wpo41UMvvE7vOOkoLKOXdB3R4P9F8em77wdmApjZVcAlwD7ApyQBTPaYI6IvITLLKvmsIm7x6r0WcPzcvfjPyWc4YMYsFozv/ra39jFxZLo0PzSz70XWrL4b3bL6tOWbsqdEGGV8KoThpjDRm9nZfZb/HvB7eeyrKMkPKvgkC8ZndBR8p30WIftBSjjdGOSKVFVQRp3ecepIXUbd1I48JV/GvoscQZRlXL3jOPWh8aIvIpuvUvLJGPKOY5hq9XU9S7aMCc6KurSgM7w0XvT9aKLkk6SJJ+S19pJ9lpOo0o6pT3vy1DCPp/epEJyyaLTo885Q6yb5FnnLPgujUL7xkTfOsNJo0fcjjfTqKnknPb3mvXGcUWSoRR9KEySfZ4x5l2+ahE9w5owijRV9v7JNXaYaroKyX3sZc984jpOdyue6qZomZPMtsoyxz5O8LiSe17w3WS8a7jhp0IwZuV+sp2wam9E7o8Mwj7ypEz4v/fDioh9S+pVvhmlMveM4vRlK0YfWqJtUtmlRRszD3iHrOKNGI0Xv2ajjOE44jRS94ziOE46L3nEcZ8hx0TuO4ww5LnrHcZySkLRM0r2S1kh6b4flH5P0o/j2U0mPJ5ZNJZatSLPfkT9hynEcpwwkjQNXAq8D1gIrJa0ws9WtNmb2R4n27wSOS2xiq5m9OMu+PaN3HMcphxOBNWZ2n5k9A1wPnNGj/dnAdXnsuJEZ/V4HTvgQS8dxSkEzZ6SZ+mORpFWJx8vNbHl8fzHwUGLZWuCkjvuUDgOOAG5OPD0n3vYk8CEz+3JoUI0UfT/2Pmh+0ElTE0fOb9xJU6Fz3fS7hmyv68c2fV4Px6mQ9Wa2tMsydXiu20UQzgJuMLOpxHOHmtk6Sc8Bbpb0YzP7eUhQXrpxHAeAzY/6r+SCWQscknh8MLCuS9uzaCvbmNm6+O99wC3sWr/viYveqQ3dZq58YnxhyZE4TiGsBJZIOkLSLCKZ7zZ6RtJRwELgu4nnFkqaHd9fBLwMWN2+bjcKE72kqyU9KuknXZZL0ifjYUZ3Sjq+qFh6UeW0v2mpOtY8pigGOk5RnAWfothpEmY2CVwI3ATcDXzBzO6SdJmk1yeang1cb2bJss7RwCpJdwDfIqrRB4u+yBr9NcDfANd2WX4asCS+nQT8LV06JjrRr0M2tE4/jAxSn8+CJjqLe3p++NzyZbHxqZlVh+CMMGZ2I3Bj23OXtD2+tMN6twEvyLrfwjJ6M7sV6FX0OwO41iJuBxZIOrCoeHpRdaYcQlkxekes4wwfVdboOw01WtypoaTzJa2StGrD09uDd9Avs01SZ9mniS3Nax5WtmyfW3UIjlMrqhR98FAjM1tuZkvNbOk+e8ze+XzeJYg6yz4vsr5nedXn68ymLZ0+ks+y8fGpnssdp65UKfo0Q40ykzbDrZPsJ46cX2o2n6Vsk7Y+360jdlhG3GzcuC247eMbns60j02PjWbfk5OdKkW/AvjdePTNrwKbzeyRtBvJO6uHesi+iBiKeK+qwkfcOE44hY26kXQdcArRKcFrgfcDMwHM7CqinufTgTXA08Bbi4olywicKs+azSL5IrP5YSnb+IgbJwsan9H1l2tTKEz0ZnZ2n+UGvCOPfYXMfdMU2Rcl+SKy+ToOq/SOWMfZnaGc66YbWWUPxV+UO2upJo9RNmVl83Wuz/friB12fPqD4WZoRB86o2XWE6mKEP6gdfhQyZeZzZdBFfV5H3HjNJmhET0UL3vYXc5pxZ9XJ2teks87m0872qYsyqjPlzHixnGyMFSih3Jkn6SKETpVZvJ5k2fZZhTq8z600slC42av1Iz8jk1NPIs0z5izZvNVdsIWUbYZ9fq8M/w0TvTQ/8SeNJlsU2S/90HzU8U6SMmmF1lq82k7YeuG1+edptNI0UP+sq+z8NPGNqjk86zNZyHvsk3d6vN1w0fcDD+NFX0IaWvUdRN+lngGrctnKdn0Is9O2GE5G9Y7Yp2yaXRn7Jz9F7Htl+t7tmmJL83FxFtyrWo++ywHm1DBFzENcZZsvg5j56FZ9XnviHWy0mjRQ5jsIXw0TpKkcIuW/iC/JPKSfFOz+SLLNl6fd2x8Ri0vopOGxoseipV9iyKkn0eZqGrJNzmbzwsfP+/UnaEQPaSTPaQr5bRThzp+mlp8FZIv6wSpQcbON6lsUxTeETsaDFVnbJr6cxNOJupGWZLPSi/JZxlSmbUTdpjKNl6fdwahcRm9ZvYOuSW2srL7MslT8NBf8mWVbKD8M2HzyuabPKzSGR0amdGHZKFps/s6Z/hp46tS8llKNr0kP0zZvNfnnapIldFL2tPMnioqmDTMXLSIHet7Z+2hdfsWSZlWneVnPfDUVfJNOQu2jhRVtvH6/OgQlNFLeqmk1cDd8eMXSfpUoZEFEJrZZxk73sqiy8z0B9ln6OscRPJZyVqyKWpIpZdtnKqQtEzSvZLWSHpvlza/LWm1pLskfS7x/Fsk/Sy+vSXNfkMz+o8Bv050nVfM7A5Jr0izo6IIyewhfXafpF28eWX7eRxE0hzEBpV8maNsqjwL1ss2ThFIGgeuBF4HrAVWSlphZqsTbZYA7wNeZmabJO0XPz9BdDnWpYABP4jX3RSy7+DSjZk9JO2SCdVmSEJLYCGlHAjrqO1FHer5eQoeipN83h2w4Nm801hOBNaY2X0Akq4HzgBWJ9q8DbiyJXAzezR+/teBb5jZxnjdbwDLgOtCdhwq+ockvRQwSbOAdxGXcepEmuy+xaDSL5u0ZagmSj5rySYP6jSkErw+30AWSVqVeLzczJbH9xcDDyWWrQVOalv/SABJ3wHGgUvN7Gtd1l0cGlSo6C8APhFveC3wdXK6sHfehMq+RV5ZftEUIXioV7lmUMrK5tPiZZtmY+Mz0nzW15vZ0i7LOn0Are3xDGAJcApwMPBtSccGrtuVINGb2XrgzaEbLRKN9w85tJSTpI5ZftYJyMqQfD+GOZsvq2zjJ0kNHWuBQxKPDwbWdWhzu5ntAH4h6V4i8a8lkn9y3VtCdxw66uazkhYkHi+UdHXAej17mCUdKulbkn4o6U5JpwfFM7EoaHTIzEWLMp352RrBUsRMj0XuO/T1hrx//SRfdsmmH57Nh+Nlm8pYCSyRdERcAj+LeIBLgi8DrwKQtIiolHMfcBNwauzehcCp8XNBhJZuXmhmj7cexL3Bx/VaIaSHGbgY+IKZ/a2kY4AbgcNDg9fEImxj/+w7S4bfol24eWb7eR1I0hzMQg6QVUi+H6OWzTvDh5lNSrqQSNDjwNVmdpeky4BVZraCZ4W+mmjAyx+b2QYASZcTHSwALmt1zIYQKvoxSQtbPcHxUJ9+64b0MBuwd3x/Prv/jOlLqOxhMOG3qCLL70bZgofiTooapGRT12x+ELxsM5yY2Y1ECW3yuUsS9w14d3xrX/dqoG8lpROhov8IcJukG+LHvwVc0WedkB7mS4GvS3onsCfw2sB4dqElsTKFXyVpy1F1kHxVJZugbRSUzXvZxqkLoZ2x10r6AVHtSMB/ayvBdCKkl/hs4Boz+4ikk4F/kHSsmU3vsiHpfOB8gEP27T6GPavwof7Sz9LXEHqW6yClGihW8oOWbJqYzTtO3qSZ6+YeYFNrHUmHmtmDPdqH9DCfRzToHzP7rqQ5wCLg0WSjeBzqcoCXLDms75CiNOWcFnWUftYphPMSPBQr+X6UUbKpYzbvZRsnb4JEH5dW3g/8kqiDQETZ+Qt7rLazhxl4mKiH+U1tbR4EXgNcI+loYA7wWK9YWpf1Gtu8oXfMKbP7JO2CLUv8g84Nn2aemjpIvsi6vLM7XrYZXUIz+j8Ajmr1/oYQ2MP8HuDvJP0R0YHj3Lgzoi8hsodd5ZdF+tBbwGkPAkVc6CNvwUO1ks8Dz+Yd51mCp0AAUn8CA3qYVwMvS7vdFi1phQgfBsvyu1GEuENJO9NkUyQ/iqNsHKdIQkV/H3CLpP8HbG89aWYfLSSqlIRm9y3yyPKrIss0wnkJHuov+VBGLZv3sk12pjSj8ddTCBX9g/FtVnyrHWmz+xZNkH7WOeLTTGEwLJLPs2TjOMNC6PDKD0C9rjDVjazCh92FWpX4B734R9mCh+IlH0LeJRvP5p1hIXTUzcnA3wN7AYdKehHw+2b29iKD64QFTGoGu8oui/Shs3DzlH/eV3PKW/BQH8nXvWTjOHUmtHTzcWp0hamWpGY/GSbdQbL8doq41N4gZJlhsk6SD6EJJRvP5p060+grTGUVPuQj/Sqpu+AhTPJl1eXTMAxTHThOkqG4wlRa4UPzpD/I3PBpLhAyrJKvawesj5t3ymCQK0yVXp/vRxbhw+4SrYP4BxF7i7wFD/WSfChpJD9M2byXbZwWoaI/ysx2ucKUpJcB38k/pMHJKvwWnSRbpPzzkHqLtJf3y1PwkI/kQ8m7Ll92B6xn805ZhIr+r4HjA54rnCmFz8M2qPCT5CnjIshy/days3jIb4RNHc5+9WzeaQo9rRkPq3wpsK+k5ET4exPNX1MJLUHtuTUsy05KMA/p14WsF+eucxZfVV2+7JKNZ/NOmfRLj2cRjZ2fAcxLPL8FOLOooEJJK3zYXY5NE39WuUO6K0C55JuLZ/P5Mq0ZA023XQd6it7M/g34N0nXmNkDAJLGgL3MbEsZAYaQRfgt6p7tDyL2FkUIHkZX8p7NO00jtOD9QUkXEI2d/wEwX9JHzewviwstPYMIHzpLtUz55yH1FlkmYcozi4f6Sz4LdZe8Z/NOJ0JFf4yZbZH0ZqJph/+USPi1En2LpOSySr9FqHx7HRDyFHg/ihQ8lN/pCsV1vA5bycZxuhEq+pmSZgJvAP7GzHZICrpASNUMmuWHUqbM28k6hWpVgodiJD/qJRvP5p1uhIr+08D9wB3ArZIOI+qQLZ1Wx8i8qU2p1muXYdHiL5JB58YuQvDgknecuhI6TfEngU8mnnpA0quKCSmMpKzSSh/yLe+UQZlyb+GSzw/P5p0q6TeO/hwz+8e2MfRJanGFqZbEsggfOku0SvnneTWbuggeqpd8Vrxk4+SFpGVE08mMA58xsw91aXcm8EXgBDNbJelwovnF7o2b3G5mF4Tut19Gv2f8d17PVjVh0Cw/ST/ZDnIgKPqyZFnH/KadUjjvLB6KlbyXbJwqkTQOXAm8jmjOsJWSVsTXzk62m0c0ceT32jbxczN7cZZ99xtH/+n47weybLxKBs3y+1HHa0jWTfDgkvds3klwIrDGzO4DkHQ9cAawuq3d5cCHgYvy2nHoFaY+DPw5sBX4GvAi4A/N7B/zCqQo2uVXlPirYNCz9bJcFGRUJJ8HfmLUcDBlY2m+K4skrUo8Xm5UMJLQAAAV2klEQVRmy+P7i4GHEsvWAiclV5Z0HHCImX1FUrvoj5D0Q6KBMBeb2bdDgwoddXOqmf2JpDfGwf0W8C2gp+hD6lGSfhu4FDDgDjN7U69ttt70+TOfDAx9V/Is71RBHqdi10nwUE/JN6Fk49l8LVlvZku7LOv0Qd85TD2edeBjwLkd2j0CHGpmGyS9BPiypOeHzlAQPI4+/ns6cJ2ZbWy72tRuhNSjJC0B3ge8zMw2SdovMJ6dssoqfOgszTrJP+/5NYoWPBSXxUOzJO/ZvNOBtcAhiccHA+sSj+cBxwK3xH49AFgh6fVmtgrYDmBmP5D0c+BIIPnroSuhol8h6R6i0s3bJe0L9PsWhdSj3gZcaWab4hfwaGA8O0nKaxDpt+gm16IOAEVPlpT1mq1FCh5c8oPi2XwjWQkskXQE8DBwFrCzgmFmm4GdZ15KugW4KB51sy+w0cymJD0HWALcF7rjvqKPf078C1HnwJZ4R08TSbsXfetRREckJH2HqLxzqZl9LTD23chb+kmaNHvdIBfkznJREJf8s7jknW6Y2aSkC4GbiHx3tZndJekyYJWZreix+iuAyyRNEs05doGZBX8Q+orezKYlfcTMTk489xTwVJ9Ve9ajEvtfApxC9DPm25KONbPHd9mQdD5wPsABBx1CCEVKv44MIneop+ChWZJ3nH6Y2Y1E84Uln7ukS9tTEvf/CfinrPsdC2z3dUm/qX6F+V3pV49qtflnM9thZr8gOhlgSfuGzGy5mS01s6ULFqYf1rh5x147b8NEHq9ry/a5tZT8xsenGid5z+aduhJao3830clTU5K2EmXrZmZ791inZz0q5svA2cA1khYRlXJ61p2mLPTY1Jl2KTYl28/7IJX1uq1ZLtpddKkGqp+J0iXv1JnQuW5SnxkbWI+6CThV0mqiutMfm1nfU05bktp79ta0Ye1GN4FWeQAo6pfHIBflLkPwUL7km1KXd5xBCD1hSsCbgSPM7HJJhwAHmtn3e63Xrx5lZkb0a6HbXDo9SYorD+knCZFt2oNBVaWjsgUPoyP5svBs3hmE0NLNp4Bp4NVEp+c+STRG/oSC4kpNkdLvRp1r/oPIHeoveKiH5L1k4zSBUNGfZGbHx6ffEp/cNKvAuAaiCulXzaBib1Gm4MEl79SfKRvL7ftVFaGi3xGf6WoA8eD96cKiypH2f9AwiT/PD19WwYNLvkg8m3fyIFT0nwS+BOwn6QrgTODiwqIqkKaKv4iMYhC5Q7MEDy55Z3QJHXXzfyT9AHgN0dDKN5jZ3YVG1oXJqXwvFN1NoFUdAMr4iViV4MEl7zhV0O8KU3OAC4DnAj8GPm1mk2UE1oukqCb23FHIPtIIt99BoQ71vUHlDtUIHkZT8p7NO3nSL6P/LLAD+DZwGnA08IdFB5WGMqTfjzqIvBN5yB0GEzy45NPiknfypp/ojzGzFwBI+nug57j5qqmD9KsmL7mDC74KXPJOEfQT/U5bxme6FhxOfrQLb1jFn6fYW1QpeKif5L0u7zSdfqJ/kaTWFUwEzI0fh8x1UyuGRfxFiL3FoIIHl/wgeDbvFEW/i4OPlxVI2XQTZl0OAEUKPUkecofqBQ8uecfpRug4+pGhn2DzOhCUJfJu1EXw4JJ3yTtF0zjRT00/K6mFe7dfx6R4qhb0IOQld6iP4MEl7xTL5JQa/b2HBoo+SVJcVUi/CeQpd8hH8FDPLB6849UZThot+iTtQhtV8ect9hZ1EjwMh+Q9m3fKYmhE386oiL8osbfIS/Dgkk/iknfKZGhF386wiL9osUO+cod8L/Pnknec9IyM6NvpJsy6HADKEHo7Lvjicck7VTCyou9GqGCzHhCqEHgv8pY75H+hbpe84wyGiz4jdRN2WkZV8OCSd6pD0jLgE8A48Bkz+1Db8guAdwBTRJdsPd/MVsfL3gecFy97l5ndFLrfxol+cvJZSU0sGNoTdwuhCLlD/oKH4ZK84wDEV+m7EngdsBZYKWlFS+QxnzOzq+L2rwc+CiyTdAxwFvB84CDgXyUdaWZBX+rGiT5JUlwu/d0pSuw7t98gwUN1kvds3ok5EVhjZvcBSLoeOAPYKXoz25Jovyfx5Vvjdteb2XbgF5LWxNv7bsiOGy36JO1SG1XxFy13KEbw4JJ3hp7FwEOJx2uBk9obSXoH8G5gFvDqxLq3t627OHTHhYq+Xz0q0e5M4IvACWa2Ko99j4r4yxD7zn254FPhkh8OktOuBLBIUtJhy81seXy/00Z2G9VhZlcCV0p6E9G1ud8Sum43ChN9YD0KSfOAdwHfKyoW6CzEJsm/TKHvtu+CBA8ueWfoWG9mS7ssWwscknh8MLCux7auB/4247q7UGRG37ceFXM58GHgogJj6UgveVZ1EKhS6EmKlDsUfwUol7xTQ1YCSyQdATxM1Ln6pmQDSUvM7Gfxw/8CtO6vAD4n6aNEnbFLSHHFvyJF37ceJek44BAz+4qkrqKXdD5wPsD8RYcWEOru1EW4ZdN0wYNL3qkn8VX6LgRuIipnX21md0m6DFhlZiuACyW9lujqfpuIyjbE7b5AlChPAu8IHXEDxYq+Z01J0hjwMeDcfhuKa1zLARY/Z2k9Tl0dIoqWe4thzuLBJe/0x8xuBG5se+6SxP0/6LHuFcAVWfZbpOj71ZTmAccCt8TXoj0AWCHp9Xl1yDrdKUvuMPxZPLjknXpTpOh71qPMbDOwqPVY0i3ARf0kPzU1zcaN25iYmFNI0MNMmXKHcgQPLnnH6Udhog+sR2UmKS2XfnfKljuMjuDBJe80g0LH0ferR7U9f0rW/bTLbJTFX4XYW5QleHDJO04ahubM2CSjJP4qxd5i1AQPLnmnWQyl6NvpJsMmHQDqIPR2yhQ8uOQdJysjIfpu9JJnFQeBOsq8nbLlDi54p1qSM+Y2lZEWfS+aIN0yqULw4JJ3nDxw0TtdqUruUB/Bg0veaT4uemc3qhQ8uOQdJ29c9A5QvdyhXoIHl7wzPLjoR5g6yB3qJ3hwyTvDReNEPzk5XXUIjaUuYm/hgneccmic6GF3YS3YZ4+KIqk/dZM71FPw4JJ3hpdGir4dF/+z1FHsLeoqeHDJO8PNUIi+nVERf52l3o5L3nGqYyhF3043ITblANAkoSeps9zBBe+MDiMh+m70E2hZB4KmirwbdRc8uOSdcFrXwGgyIy36fgybgIumCYIHl7wzerjonYFxwTtOvXHRO5lpiuDBJe+MNi56JxVNknsLl7wz6rjonSBc8I7TXFz0TleaKPcWLnnHeZaxqgNw6sWmxzbvvDWRzY9udMk7tUXSMkn3Sloj6b0dlr9C0n9ImpR0ZtuyKUk/im8r0uzXM3qnsVJvxwXv1BlJ48CVwOuAtcBKSSvMbHWi2YPAucBFHTax1cxenGXfhYpe0jLgE8A48Bkz+1Db8ncDvwdMAo8B/8PMHigyJmd4xN7CBe80hBOBNWZ2H4Ck64EzgJ2iN7P742W5TtNbWOkmcfQ6DTgGOFvSMW3NfggsNbMXAjcAH+633anJqaETVRk0vSTTDZe8UzMWSVqVuJ2fWLYYeCjxeG38XChz4m3eLukNaYIqMqMPOXp9K9H+duCc0I23C2vhvvMHiXWoGDaZd8IF75TF5OR0mrPk15vZ0i7L1OE5SxHKoWa2TtJzgJsl/djMfh6yYpGi73T0OqlH+/OAr3ZaEB8VzweYO++gjit3ktuoyH8UxN7CBe80mLXAIYnHBwPrQlc2s3Xx3/sk3QIcB1Qu+uCjl6RzgKXAKzstN7PlwHKAhfu/IPgIOGzyHyWhd8Il7zSclcASSUcADwNnAW8KWVHSQuBpM9suaRHwMgJK3S2KFH3Q0UvSa4E/A15pZtsLjAfoL8uqDwSjLvNOuOCdYcDMJiVdCNxENEDlajO7S9JlwCozWyHpBOBLwELgv0r6gJk9Hzga+HTcSTsGfKhttE5PihR936OXpOOATwPLzOzRAmMJxkVbH1zwzrBhZjcCN7Y9d0ni/kqipLh9vduAF2Tdb2GiDzl6AX8J7AV8URLAg2b2+qJicpqBC95x8qXQcfQBR6/XFrl/p1m44B2nGPzMWKdyXPCOUywueqcyXPCOUw4ueqd0XPCOUy4ueqc0XPCOUw0ueqdwXPBOkxmG+bVc9E5huOAdpx646J1ccbk7Tv1w0Tu54IJ3nPrioncy43J3nGbgondS44J3nGbhoneCcLk7TnNx0Ttdcbk7znDgond24mJ3nOHERV8A27du4OktD7PH3ouZPXefqsPpicvdcYYfF33OPHTvv/DDm/8XGpuBTU9y3Gs+yCFH/kbVYe3Exe446ZjaMdn4703jRN960+fvN1F1KLuxfesGfnjz/2JqctvO5374zfex3yEnV5bZN/0D6jjO4DRO9C06Caxq+T+95WE0tutbqrEZPL3l4VJE71J3HKcTjRV9J7qJrqwDwPjMPZieemaX56andvDMts1s37ohV9m71B3HCWWoRN+NXlLsdRBI06naqs0bAmBsxhxsegqbnuL7X3tX6nq9i9xxnLwYCdH3optQ//OBm7hn5RU7O1Wfd8LFHHDYqcDuB4dOtXmbngKE2Q4mn3kSgP/41/cyZ87RzJqzsJgX4ziO04GRF30nntm2iXtWXsH01HaY2g7A3d+/nHkLlrDn/CN2OTg8s20TG9Z9BxjbZRvSOCBsl2fHeGLTvcxbeBTbnnqEOXse6NJ3HKdwRlr0z2zbxLanHmF8xlymJrfuFO+2px6JOlVjyQPY9DN8/+u/w9EnXsLE/iew7alH2LLpHtb86OOgcaYnn95l22bTSNrluemprdxx67vR2Bhj47N2+6XgOI5TBIWKXtIy4BPAOPAZM/tQ2/LZwLXAS4ANwH83s/uLjKlFqzRjZtj0M4yNzwbgeSdczMT+J2DTk7utY9M7WP29Sxkbm9FR7gDjM/bAbIrnnXAxAPes/PPol8FOorr91PSOncsn9j/BM3vHGQEGcaKk9wHnAVPAu8zsptD9jvVvkg1FtYsrgdOAY4CzJR3T1uw8YJOZPRf4GPAXRcWTJFmaselolMz01Hamp7Zzz8o/ByLha2zW7ivbVNS2g+THZuzBkce9h5f+xj9zwGGncsBhp/KCl32YsfE5XWPR2Ay2PfVIPi/McZzaMogT43ZnAc8HlgGfircXRGGiB04E1pjZfWb2DHA9cEZbmzOAz8b3bwBeo/Z6RwHsLM10oCXeAw47lRNPvRaNzQzfsE2xz0Ev2yU7n7fwKGir1O+yyvQkc/Y8MHwfjuM0lUGceAZwvZltN7NfAGvi7QVRZOlmMfBQ4vFa4KRubcxsUtJmYB9gfbKRpPOB8+OH22/+/Ek/GTC2GcAL6HCgm9rx1PSqf33rj4FW7WYhcDiRrZW4JZkGND21/Rf//s/LNnXYX3IbY23r3N9lnTxZRNt7WhM8rvTUNba6xnXUoBt4YtM9N938+ZMWBTafI2lV4vFyM1se3x/EiYuB29vWXRwYU6Gi75SZt6e2IW2I36jlAJJWmdnSwcPLl7rGBfWNzeNKT11jq3Ncg27DzJblEQuDOTHIld0osnSzFjgk8fhgYF23NpJmAPMBP1PIcZxhZBAnhqzblSJFvxJYIukISbOIOhJWtLVZAbwlvn8mcLOZBR+lHMdxGsQgTlwBnCVptqQjgCXA90N3XFjpJq4vXQjcRDSU6Gozu0vSZcAqM1sB/D3wD5LWEB21zgrY9PL+TSqhrnFBfWPzuNJT19g8rj4M4sS43ReA1UT9h+8ws6nQfcsTaMdxnOGmyNKN4ziOUwNc9I7jOENObUUvaZmkeyWtkfTeDstnS/p8vPx7kg6vSVzvlrRa0p2SvinpsDrElWh3piSTVNpQuJDYJP12/L7dJelzdYhL0qGSviXph/H/8/SS4rpa0qOSOp4voohPxnHfKen4msT15jieOyXdJulFdYgr0e4ESVOSziwjrlphZrW7EXVU/Bx4DjALuAM4pq3N24Gr4vtnAZ+vSVyvAvaI7//PusQVt5sH3Ep04sXSGv0vlwA/BBbGj/erSVzLgf8Z3z8GuL+k9+wVwPHAT7osPx34KtHY6l8FvleTuF6a+B+eVpe4Ev/vm4EbgTPLiKtOt7pm9HWdPqFvXGb2LTNrTYRzO9F416IJeb8ALgc+DGzrsKzK2N4GXGlmmwDM7NGaxGXA3vH9+aQYtzwIZnYrvc8nOQO41iJuBxZIKnwejX5xmdltrf8h5X32Q94vgHcC/wSU8dmqHXUVfadThdtP993lVGGgdapw1XElOY8o8yqavnFJOg44xMy+UkI8SULesyOBIyV9R9Lt8Qx/dYjrUuAcSWuJMsF3lhBXCGk/h1VQ1me/L5IWA28Erqo6lqqo63z0uU2fkDPB+5R0DrAUeGWhEcW76/DczrgkjRHNhHduCbG0E/KezSAq35xClAV+W9KxZvZ4xXGdDVxjZh+RdDLR+OZjzWy6wLhCqOKzH4ykVxGJ/uVVxxLzceBPzWyqhDkTa0ldRZ/mVOG1JU6fEHQasqTXAn8GvNLMtrcvryCuecCxwC3xB/0AYIWk15vZwHOBDBhbq83tZrYD+IWke4nEv7LiuM4jmhIWM/uupDlEk3dV/fN/oNPhi0TSC4HPAKeZ2Yaq44lZClwff/YXAadLmjSzL1cbVolU3UnQ6UZ0ALoPOIJnO8qe39bmHezaGfuFmsR1HFEn35I6vV9t7W+hvM7YkPdsGfDZ+P4iorLEPjWI66vAufH9o4lkqpLet8Pp3un5X9i1M/b7JX7WesV1KNH0uS8tK56QuNraXcMIdsbWMqO34qZPKCOuvwT2Ar4YZxAPmtnraxBXJQTGdhNwqqTVRFfP+WMrOBsMjOs9wN9J+iOi0si5FtuiSCRdR1TGWhT3D7wfmBnHfRVRf8HpRFJ9Gnhr0TEFxnUJUT/Zp+LP/qSVMKNlQFwjj0+B4DiOM+TUddSN4ziOkxMuesdxnCHHRe84jjPkuOgdx3GGHBe94zjOkOOid2pLPNPgjyTdIek/JL00fv4gSTd0WedwSW9KPD5X0t+UFbPj1BEXvVNntprZi83sRcD7gA8CmNk6M9ttqtn4DOnDgTe1L3OcUaaWJ0w5Tgf2BjZBlLUDXzGzYyWdS3Sm6BxgT2AP4GhJPyKa3XQTcJCkrwG/AnzJzP6k9Ogdp0Jc9E6dmRsLew5wIPDqLu1OBl5oZhslnQJcZGa/AVHpBngx0dQU24F7Jf21mT3UZVuOM3R46capM63SzfOI5sO5tss1B75hZr0mtPummW02s23AaqCUq345Tl1w0TuNwMy+SzTh2b4dFj/VZ/XkDKJT+C9ZZ8Rw0TuNQNLziCYf6zfZ2RNE0zI7jhPjmY1TZ1o1eoim5H2L9b94xJ3ApKQ7iKak3dSrseOMAj57peM4zpDjpRvHcZwhx0XvOI4z5LjoHcdxhhwXveM4zpDjonccxxlyXPSO4zhDjovecRxnyPn/udfgHLaexu0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create grid for computing posterior intensity\n",
    "b = np.linspace(0,1.4,50)\n",
    "p = np.linspace(0,1.4,50) \n",
    "B,P = np.meshgrid(b,p)\n",
    "\n",
    "# evaluate posterior intensity over grid\n",
    "I = list()\n",
    "for ind in range(len(P)):\n",
    "    l = list()\n",
    "    for i in range(len(P)):\n",
    "        l.append(posterior.eval([B[ind][i],P[ind][i]]))\n",
    "    I.append(l)\n",
    "I = I / np.max(I)\n",
    "\n",
    "# plot posterior intensity\n",
    "plt.figure()        \n",
    "plt.style.use('seaborn-bright')\n",
    "plt.contourf(B,P,I, 20, cmap = 'coolwarm', vmin = 0, vmax = 1)\n",
    "plt.title('Posterior Intensity')\n",
    "cb = plt.colorbar()\n",
    "cb.ax.set_title('Scal. Intens.')\n",
    "plt.gca().set_ylabel('Persistence')\n",
    "plt.gca().set_xlabel('Birth')\n",
    "plt.xlim([0,1.4])\n",
    "plt.ylim([0,1.4])\n",
    "plt.scatter(pd_t[:,0],pd_t[:,1], s = 20 ,color = 'black') # overlay persistence diagram\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine Learning Example: Classification\n",
    "\n",
    "Next, we'll give an example of how our Bayesian framework can be used in machine learning problems. We'll consider the problem of distinguishing circles from ovals. First, we define some auxillary functions to create our datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# draw samples from circle\n",
    "def noisy_circ(n_pts,noise_lev = 0.1):\n",
    "    t = np.linspace(0,1,n_pts)\n",
    "    x,y = np.cos(2*np.pi*t), np.sin(2*np.pi*t)\n",
    "    circ = np.array([[h,v] for h,v in zip(x,y)])\n",
    "    noise = np.random.normal(scale = noise_lev,size = [n_pts,2])\n",
    "    n_circ = circ+noise\n",
    "    return n_circ\n",
    "\n",
    "# draw samples from oval\n",
    "def noisy_oval(n_pts,noise_lev = 0.1):\n",
    "    t = np.linspace(0,1,n_pts)\n",
    "    x,y = 1.5*np.cos(2*np.pi*t), np.sin(2*np.pi*t)\n",
    "    oval = np.array([[h,v] for h,v in zip(x,y)])\n",
    "    noise = np.random.normal(scale = noise_lev,size = [n_pts,2])\n",
    "    n_oval = oval+noise\n",
    "    return n_oval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll create 100 noisy circle and ovals, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAC7CAYAAABxX7W7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAE9ZJREFUeJzt3X9oHdeVB/Dvaaxu46S2BXaJiGy3sFZcVywKeql3cchCV0u0xd2mhUIbaFLWRSzEoEIJbQn+Z/1noCC2gcXEIQ007S50101S1m4FBdOlK6yXfaxlO1GyBUVytXRNnyRT/6iUnv1DGmX0NPNm5s19c+/c+X4gJPqReVfzzpx35sy9M6KqICIif3zI9gCIiMgsJnYiIs8wsRMReYaJnYjIM0zsRESeYWInIvJM7sQuIvtF5Bcick1ErojIuImBEdnG2Kaykrzz2EWkD0Cfqr4pIh8FUAfwhKpeNTFAIlsY21RWuSt2VV1U1Tc3/vsmgGsAHsy7XSLbGNtUVkZ77CLycQAPA5gyuV0i2xjbVCY7TG1IRO4H8GMA31DVlYifjwEYA4D77rtv+PDhw6ZemmiLer1+Q1X3mdoeY5tckTa2c/fYAUBEegC8AeCCqn436fdrtZpOT0/nfl2iKCJSV9WaoW0xtskZaWPbxKwYAXAWwLU0gU9UFoxtKisTPfZjAL4K4DMi0tj457MGtktkG2ObSil3j11VfwlADIyFyCmMbSorrjwlIvIMEzsRkWeY2ImIPMPETkTkGSZ2IiLPMLETEXmGiZ2IyDNM7EREnmFiJyLyDBM7EZFnmNiJiDzDxE5E5BkmdiIizzCxExF5homdiMgzTOxERJ5hYici8gwTOxGRZ5jYiYg8w8ROROQZJnYiIs8wsRMReYaJnYjIM0zsRESeYWInIvIMEzsRkWeY2InIqvpcE0+dnUJ9rml7KN5gYiciqyYmZ3HxnRuYmJy1PRRvMLE7LlzNsLIhH42PDOCxQ3sxPjIAYHsFz7jPboeJjYjISwCOA/itqg6a2GYV1OeamJicxfjIAIYP9kb+TlDNBIL/fuXE0ULGWGWM62IMH+zdEs/hmH/lxNFtX1MyUxX7ywBGDW2rMtKcggbVzOhgH1burGGof3dsZUPGvQzGdWrt4jFLrLZW8K1fUzIjiV1VLwL4nYltVUm7gA0OBGC9Sjk/s4jG/BJ23duzWd2zN9ldjOv06nNNfP37lzbjsTWRZ4nVoIIP4nz4YC/GRwY2t5tnjFUphNhjT9DNYGgN4LDWAyHqQ4CVDJnWabxPTM6ieWsVvTt7NpNwUvxm3X7eIqZKhZCRHnsaIjIGYAwADhw4UNTL5marvxccAKODfXjq7BTGRwa2vX5rbzJJmp4+ZVfW2I5y+o2raMwvYeXOGs49cyz1/xdumwQVdvj7WWO13fazxHH4d1vH5LPCKnZVPaOqNVWt7du3r6iXzS1PpRFX/aS56h8cCOdnFo1VGVWqWIpU1tiOpLr13ylFtU/izkY7Ed5eljgO/67pMbmssIq9rPJUGnHVfvD9y9eX8eLTj7Q9KzBZZVSpYqHOnPrcpzYrXFdlieOqxryRil1EfgjgVwAeEpEFETlhYrtlNzrYh96dPRgd7Nvy/fGRAfTu7EHz1urmQRR3VtCuyuikH7pyZw2nX79SiQtIeVUxrstQ1WYZYxn+nm4wNSvmK6rap6o9qtqvqmdNbLfszs8sonlrFednFgFsneny4tOPbCbzToMva2tlYnIWjfklNBaW2Y5JgXFNZcVZMRmk6Y2Hv9daiZvo97XbfpLxkQEM7d+zZS48UZjPK519+3vaYY89gzQr4lq/F+6Zm+j3tdt+kuGDvZlmOlD1+LzSuUorWJnYM2hNzFGJul3yzjvlq932OZWRTGgX02VXpQupohmnNZlQq9V0enq68Nf12VNnp3DxnRt47NBe76uRJCJSV9WajdeuWmyXvaAo2/jTxjZ77J7gKlSywZW1EXlWzLowftOY2GOU7UJLVad1UTam47pbBUXWcSYl6Ljt+VoQscceowoXWsp2Gkr5mY5rE9eNomQdZ1L/PG573Rq/bazYY/j6SR7WSZVTtjMZWhe8b6ODfaWI66zHX9IZa9L28sa1a8cFK3ZEV66+fpKHJc2wWbmzhsb8EoAPqpzTr19BY2EZK7dXce7ko8UOmFKJiueynYGaPv6SthfeP8HdKbOcybq2f5nY4d6bUpS4YA/2x1D/7m2PLPuf//v9+i+JFDlUyiAqnqs01a8T4f3TST5wbf8ysSP6Taly/7n1FqyBiclZ3Ly7ht6dPTh1/Iit4RHax+foYB8uX1/eco+iKpyB5hHeP64l6U4wsSM66KtaxQPxSSAu4VNxNttkt1fRWFgGsD0+w/coevJoue8Pb0MnH4Ku5Qsm9hjtPrWrWs2z6rMvSCCHPnZ/5J1DAT8qTpdFHf+u7XPOiomR5bF1VeXaTIAqGB8ZwFD/bvzv8p0tdw4N45qG7oo6/l3b56zYM6rPNbFyexVD+/c48+lsi2unn1UwfLAXu+7t2bzWUfUYtMG16jwKE3tGE5OzaCws47FDe535dLalDAHuI17rsKsMLUkm9oyqmMzirimUIcB9xP1OSdhjz8i1XloRwj1F9tXdwfeC4jCxU6LwcuzWC0dMLvbwIr4daZ6kZhtbMRGqOp0xTrvFG7yAak8V24IuSPMkNduY2CO4+Ea5orW/y+RiD3vtdrTGfNRKX9vYiolQhTs7mlLFaw5FcfEUn7bHfHilb6A+18QT3/slnnjhP6y8f0zsEcJvXFUPrjL0EX3HHno5RBWCwbToxvySlfePrZgEVW3LtP7dvF1vMcLXd9jmKoeoltj4yABWbq8CIlbeP1bsIUFV+urUe5vV6ehgX+w9OXy2rQoJbtMbul0vq3jzwlV6OGFwP5fL8MFenDv5KM49c8zKmT8r9pDgoLp8fRnNW6ub36/infJaq5BTx4/g9OtXAKwn9OGDvZU9m+mmqCqd+7n88j7IIysm9pDwVe7zM4tbDq6qnw4H9ygJqslXThxlq6AL4k7rw/+m8sn7II+sRFW7suF2arWaTk9PF/66UThnPb2y7CsRqatqzcZrm4jtsuxn6kye9zdtbFe+x86ZB+lxamMxuLrXb0UcR0YSu4iMisjbIvKuiHzbxDaL0u2nl1O52Yjt1phk8UFZ5e6xi8g9AF4A8NcAFgBcEpHXVPVq3m0XIcvTy3nhqlpsxTZX9/qhPtdcn3AgglPHjxR6pmvi4umnAbyrqr8GABH5EYDPAyhFYm/V2v/iQVVpTsQ2bx1QTsEiJQA4/cZV7PrIjsKum5hI7A8CmA99vQCgdFG4+ZDgO2tozC8BwGYfjAdVZXkR22RHeJESVAs98zeR2CXie9um2ojIGIAxADhwwL354EHLZah/N+8TQwEvYpvs2XVvz5ZrJUXlFROJfQHA/tDX/QB+0/pLqnoGwBlgfUqYgdc1io8bSydpqpZnU/W8iG0qVtzZf5Fn/iZmxVwCcEhEPiEiHwbwZQCvGdhuoTiVL52kqXiezeBwJrY5O6s8NidcqFo7+89dsavqmoicBHABwD0AXlLVK7lHRk5KetCGTxebXYptzs4qDxfO/iu98jSqbeBZK6HrXNxfZV95GgjvWwDO7WcqXtrYrvS9YqKqIFZG2XDWUPe0xiL3M6VV6cQe1TbwqZVA5cZYpE5VuhVDfvKlFRPHxfYXFYM3ATOEsxGoCFnizLOZR16ynTeY2BPwIGrPdgD7Ikuc8WHr7rOdNyrdY09jdLAPl68vV+7ReGnxYrMZ4X56UquFF6zdZ/v6CCv2BOdnFjcfjUfrwlU6q0czwgvkbFd7lJ/t59WyYk9g+5PXRZyG112MufKJO8uydUbLij0kql/MWw1sNz4ygKH9e7Bye5W99YLxmoab4s6ybJ3RMrGH8BQ4neGDvdj1kR1oLCxzX3VBuzhkjLopLoHbKgzZignhKXB63Ffd027fcr+7ybUL2lygtIGLPtor0/4p2wKlMu1bai94L0cH+3B+ZtH4e8p7xWTEaXvtcf90T6f7lh8I7mi9B/vl68to3loFYOd4YY99A6fttRfsn9HBPl68Myxt7Hl+73unJV20br0H+7OPH7aaT1ixb0jqkVW5Ogr/7azczUvbn/X53veuS4r7qHuwP3nU3mMSmdjbCPfLnr/wltVTK5vCQc1kYk/rvnftgp3PkuLetfeCrZg2goQWJPXenT1eJrSo08y41aWc128P931xWo+J1n3f+nPX1hewYm8jSOLdusJtW9xDdwHg9BtX0ZhfwsqdNZx75phT1Qi11+2ZGVWQ1Hpp/blrLUom9jbCp1c2+2XdEgTjUP/u7Rd6gmmwFqbDUufqc0187aUp3Lz7Pv7rvSXcvLsGwI1kUyZJrZfWn7vWomRir6j6XBMrt1cxtH8PTh0/sq2iO/W5T2153ia5J+qC/sTkLG7efR8A8MCuP8HDe/bwPexAUs+89eeu9diZ2CtqYnIWjYVlPHZoL28NW1JRp//jIwNYub0KiER+YFM1MLHnUJ9r4vTrV2IPIpenSLp26kjZRb2Hwwd7ce7ko5tfuxyD1D2cFZNB1AKRxsIyGvNLpbthU9QMC9eu7FN7aWbJdBKDvsSBL39HJ1ixZxC1QGTl9ip+/4f3sXJnDfW55paDrAxVMRcf+a2TGPQlDnz5OzrBxJ5B1AKRcycfxVNnpzarItcuqCSdinPxkd86iUFf4iDr3+HTNFEm9gyiDpL6XBMrd9Yw1L/byQMh61Jo2x9EZJ8vcZD1NiHBsWL7Bl4mMLHnNDE5i8b8UuzsEhvCAVu2pdBERajPNfH171/aksCjFiSWFRN7Ti6etrZ7JilnSVBePsTQxOTsttuE+LQgkbNicnLx/h3tbgPr8kwdsivtLBITMWR7xkpwjLz49CNOHbum5ErsIvIlEbkiIn8UEStPrHGR7aBt92HD+86nU8XYTpuwW2Ook3i3XWC4WJCZlLdinwHwRQAXDYzFG3FBG3cApD0wkn4vzXZ8D2iDKhfbUR/6UTHVGkOdJOksBYbtQqmMciV2Vb2mqm+bGowv4oI27gBIe2Ak/Z7tKsgnVYztqA/9NDEVjvd2SfjVqffw8D/8DK9OvZepwEgb1/wA+AAvnnZB3EyTuAutrd+PuziV5o5zK7dXIxdLEXUizeSAcLwHazqA7VMFg+caPH/hrbYXJ1vjP+0EhSovSGqVmNhFZBLAAxE/ek5Vf5L2hURkDMAYABw4UO4rzp2KS/it3w/Ppw1f3ImbRx8+CHbd2xO5WIq2Y2wnyzodtl0Sfvbxw3j+wlt49vHDbbfRmqDTjmF0sA+Xry9jdLAv9Xh9lZjYVXXExAup6hkAZwCgVqt5eZNvU9PAjvTtwsV3bqB5azUxQfM5mJ1jbJsTjv24eH3y6IFU0wg7jeHzM4to3lrF+ZnF0k9XzIutGINMnQr+8/Q8AGDHhyQxuJOeg+nDnGPKr9txYLIN0umiORY1H8iV2EXkCwD+EcA+AD8VkYaqPm5kZCVkKrDCp6xJB2HSQcC+Y2d8i+1ux4ELSZWrqD8gauHRZ7VaTaenpwt/3SqqYsUuInVVtTL33NXYrmIc+ChtbLMVUxKdHpisYghgHFQNbylQEpyjTkRpsWIvCRd6mEQuSXo0ZZWxYi8J3gqAaKukR1NWGSt2IiqlYKU1JHlacNUwsRNRKQWPpqTt2IohIvIMEzsRkWeY2ImIPMPETkTkGSZ2IiLPMLETEXmGiZ2IyDNM7EREnmFiJyLyDBM7EZFnmNiJiDzDxE5E5BkmdiIizzCxExF5homdiMgzTOxERJ5hYici8gwTOxGRZ5jYiYg8w8ROROQZJnYiIs8wsRMReYaJnYjIM7kSu4g8LyJvich/i8i/icgeUwMjsomxTWWWt2L/OYBBVf0zALMAvpN/SEROYGxTaeVK7Kr6M1Vd2/jyPwH05x8SkX2MbSozkz32vwPw7wa3R+QKxjaVyo6kXxCRSQAPRPzoOVX9ycbvPAdgDcAP2mxnDMDYxpd3RWQm+3C7bi+AG7YHEYHjyuahNL/E2HYCx5VNuthW1VyvIiJPA/h7AH+lqrdS/j/TqlrL9cJdwHFl4/u4GNvdx3Flk3ZciRV7wouMAvgWgL9MG/hEZcDYpjLL22P/HoCPAvi5iDRE5J8MjInIBYxtKq1cFbuq/mmH/+uZPK/bRRxXNt6Oi7FdGI4rm1Tjyt1jJyIit/CWAkREnrGW2F1dsi0iXxKRKyLyRxGxflVcREZF5G0ReVdEvm17PAAgIi+JyG9dm9YnIvtF5Bcicm3jPRy3MAbGdbrxOBfXgJux3Ulc26zYXV2yPQPgiwAu2h6IiNwD4AUAfwPgCICviMgRu6MCALwMYNT2ICKsAfimqn4SwJ8DeMbC/mJcJ3A4rgE3YztzXFtL7K4u2VbVa6r6tu1xbPg0gHdV9deq+gcAPwLwectjgqpeBPA72+NopaqLqvrmxn/fBHANwIMFj4FxnczJuAbcjO1O4tqVHjuXbEd7EMB86OsFFJyoykpEPg7gYQBTFofBuI7GuO5Q2rjONd0xxSCMLNm2MS5HSMT3OI0pgYjcD+DHAL6hqitd2D7jOh/GdQeyxHVXE7uqjrT7+caS7eNYX7Jd2BubNC6HLADYH/q6H8BvLI2lFESkB+vB/wNV/dduvAbjOjfGdUZZ49rmrJhgyfbfcsl2rEsADonIJ0TkwwC+DOA1y2NylogIgLMArqnqdy2NgXGdjHGdQSdxbbPH7uSSbRH5gogsAPgLAD8VkQu2xrJxEe4kgAtYv2DyL6p6xdZ4AiLyQwC/AvCQiCyIyAnbY9pwDMBXAXxmI6YaIvLZgsfAuE7galwDzsZ25rjmylMiIs+4MiuGiIgMYWInIvIMEzsRkWeY2ImIPMPETkTkGSZ2IiLPMLETEXmGiZ2IyDP/Dy/LyedX6JQhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_pts, data_sz = 100,100\n",
    "circs,ovals = [noisy_circ(n_pts) for i in range(data_sz)], [noisy_oval(n_pts) for i in range(data_sz)]\n",
    "\n",
    "# visualize the data\n",
    "c,o = circs[0], ovals[0]\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "ax1.scatter(c[:,0], c[:,1], s = 2)\n",
    "ax1.set_ylim([-2,2])\n",
    "ax1.set_xlim([-2,2])\n",
    "ax1.set_aspect('equal')\n",
    "ax2.scatter(o[:,0], o[:,1], s = 2)\n",
    "ax2.set_ylim([-2,2])\n",
    "ax2.set_xlim([-2,2])\n",
    "ax2.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we compute the Rips filtrations and 1-dimensional persistence diagrams (in birth-persistence cooridnates) for each of these point clouds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pd_tilted(point_cloud, hom_dim = 1):\n",
    "    st = gd.RipsComplex(points=point_cloud, max_edge_length=3.).create_simplex_tree(max_dimension=hom_dim + 1)\n",
    "    dgm = st.persistence()\n",
    "    pairs = st.persistence_pairs()\n",
    "    simps = [s for s in filter(lambda x: len(x[0]) == hom_dim + 1,pairs)]\n",
    "    pd = np.array([[st.filtration(x),st.filtration(y)] for x,y in simps])\n",
    "    pd_t = np.array([[b,d-b] for b,d in pd])\n",
    "    return pd_t\n",
    "\n",
    "circ_diags, oval_diags = [pd_tilted(pc) for pc in circs],[pd_tilted(pc) for pc in ovals]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a persistence diagram $\\mathcal{D}$, we estimate the Bayes factor for the competing hypotheses '$\\mathcal{D}$ corresponds to a noisy circle sample' and '$\\mathcal{D}$ corresponds to a noisy oval sample', i.e. $BF(\\mathcal{D}) := \\frac{\\mathbb{P}(\\mathcal{D}|\\mathrm{Circle})}{\\mathbb{P}(\\mathcal{D}|\\mathrm{Oval})}$. Large and small values of $BF(\\mathcal{D})$ provide strong evidence for the former and latter hypotheses, respectively. Although we can use our Bayesian model to compute the numerator and denominator of $BF(\\mathcal{D})$, for numerical reasons, we choose instead to work with $\\log BF(\\mathcal{D})$. We estimate $\\log \\mathbb{P}(\\mathcal{D}|\\mathrm{Circle})$ by computing $\\log \\mathbb{P}(\\mathcal{D}|\\mathcal{T})$ using our Bayesian model, where elements of $\\mathcal{T}$ are persistence diagrams associated with noisy circle samples. The Bayes factor can be used for classification, as summarized in the following schematic.\n",
    "<img src=\"gudhi_fc.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function computes the log posterior probability $\\log \\mathbb{P}(\\mathcal{D}|\\mathcal{T})$. The function also requires the hyperparamters for computing posteriors in our framework: a prior intensity, an intensity for unanticipated features, and a covariance matrix for the likelihood. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "def training_posterior(D,T,prior,unant,lik_cm):\n",
    "    posts = [Posterior(prior = prior, clutter = unant, Dy = diag, sy = lik_cm) for diag in T]\n",
    "    pp_log = 0\n",
    "    for feat in D:\n",
    "        post_pri = [p.eval(feat) for p in posts]\n",
    "        pp_log += np.log10(np.sum(np.array(post_pri)))\n",
    "    return pp_log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use a flat prior (i.e., one with a high covariance magnitude) and drop the unanticipated intensity by setting its weight to 0. The high uncertainty in the prior and lack of an unanticipated component lead to posteriors that are entirely data-driven. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_pri = np.array([[0.2,1]]) # prior mean\n",
    "w_pri = np.array([1]) # prior weight\n",
    "sig_pri = np.array([20]) # prior covariance magnitude\n",
    "\n",
    "pri_train = Prior(weights = w_pri,mus = mu_pri, sigmas = sig_pri)\n",
    "unant = Prior(weights = np.array([0]),mus = mu_pri, sigmas = sig_pri) # unanticipated intensity set to zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create our validation set by drawing 20 more samples each from a circle and an oval then computing persistence diagrams."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_pts, data_sz = 100,20\n",
    "circs_val,ovals_val = [noisy_circ(n_pts) for i in range(data_sz)], [noisy_oval(n_pts) for i in range(data_sz)]\n",
    "c_dgms_val, o_dgms_val = [pd_tilted(pc) for pc in circs_val],[pd_tilted(pc) for pc in ovals_val]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute the log Bayes factors for each diagram in our validation set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "lbf_circs,lbf_ovals = [training_posterior(diag,T=circ_diags,prior = pri_train,unant = unant,lik_cm = 0.1)- \\\n",
    "                     training_posterior(diag,T=oval_diags,prior = pri_train,unant = unant,lik_cm = 0.1) for diag in c_dgms_val],\\\n",
    "                    [training_posterior(diag,T=circ_diags,prior = pri_train,unant = unant,lik_cm = 0.1)- \\\n",
    "                     training_posterior(diag,T=oval_diags,prior = pri_train,unant = unant,lik_cm = 0.1) for diag in o_dgms_val]\n",
    "                     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we visualize the Bayes factors to get a sense of how well our model distinguished the two classes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHYZJREFUeJzt3XmUFfWZ//H3RxZbFDQiMUbsNGaihk0ITY5BjwqJxAWMGol6NKMmGbI4BtQxIRoTlZnoqDGGxN9vgg4xZnEjccE1G4uOqIgggkuMW+wRI6IgqCjLM39UtWmh+3b17arb3dzP65x7um7dWp4vDQ91v/Wt56uIwMzMtn7bdHQAZmZWGU74ZmZVwgnfzKxKOOGbmVUJJ3wzsyrhhG9mViWc8M3MqoQTvplZlXDCNzOrEt07OoCmdtlll6irq+voMMzMuoyFCxe+GhH9smzbqRJ+XV0dDz/8cEeHYWbWZUh6Ieu27tIxM6sSTvhmZlXCCd/MrEp0qj58M7NG69evp6GhgXXr1nV0KJ1CTU0N/fv3p0ePHmUfwwnfzDqlhoYGevfuTV1dHZI6OpwOFRGsXLmShoYGBgwYUPZxCu3SkbSTpJmSnpT0hKRPFXk+M9t6rFu3jr59+1Z9sgeQRN++fdv9bafoK/wfA3dHxLGSegK9Cj6fmW1FnOz/IY8/i8ISvqQ+wIHAKQAR8S7wblHnMzOz0oq8wt8TWAH8XNK+wEJgUkS8WeA5zWwrNev+fI83flTr27z88stMnjyZBQsWsO2221JXV8eUKVOYNm0aM2fOzHyuU045hXHjxnHssce2I+L2KzLhdwc+AZweEQ9K+jEwBTiv6UaSJgITAWprawsMx6xzyTuBNcqSyKx1EcHRRx/NySefzPXXXw/A4sWLWbNmTbPJfsOGDXTv3rnHwRR507YBaIiIB9P3M0n+A3ifiJgeEfURUd+vX6ZyEGZmhZs9ezY9evTga1/72nvrhg0bxh577MHgwYMBuOaaa5gwYQLjx49n7NixAFxyySUMGTKEfffdlylTpmxx3IULF3LQQQcxYsQIPvvZz7J8+XIApk2bxsCBAxk6dCjHH398IW0q7L+jiHhZ0ouS9o6Ip4BPA48XdT4zszwtXbqUESNGtLrd/PnzWbJkCTvvvDN33XUXt9xyCw8++CC9evXitddee9+269ev5/TTT+fWW2+lX79+3HDDDZx77rnMmDGDiy++mOeee45tt92WVatWFdKmor9/nA78Oh2h8yxwasHnMzOrqEMOOYSdd94ZgD/+8Y+ceuqp9OqVDEhsXN/oqaeeYunSpRxyyCEAbNy4kd122w2AoUOHcuKJJ3LUUUdx1FFHFRJroQk/IhYD9UWew8ysCIMGDcp0Y3b77bd/bzkiSg6fjAgGDRrE/Pnzt/jsjjvuYN68edx2221MnTqVZcuW5X5PwLV0zMyaMWbMGN555x2uuuqq99YtWLCAF15ouRrx2LFjmTFjBm+99RbAFl06e++9NytWrHgv4a9fv55ly5axadMmXnzxRUaPHs0ll1zCqlWrWLt2be5t6ty3lM3MUpUefSSJm2++mcmTJ3PxxRdTU1NDXV0dV1xxRYv7HHrooSxevJj6+np69uzJ4Ycfzg9+8IP3Pu/ZsyczZ87km9/8JqtXr2bDhg1MnjyZvfbai5NOOonVq1cTEZxxxhnstNNO+bcpInI/aLnq6+vDE6BYtfCwzNKeeOIJPv7xj3d0GJ1Kc38mkhZGRKauc3fpmJlVCSd8M7Mq4YRvZlYlnPDNzKqEE76ZWZVwwjczqxIeh29mXcKsd+/L9Xjjex7Q6jYNDQ2cdtppPP7442zatIlx48Zx6aWX0rNnzzada86cOVx22WXcfvvt5YabC1/hm5k1IyI45phjOOqoo3j66af5y1/+wtq1azn33HM7OrSyOeGbmTXjz3/+MzU1NZx6alLzsVu3bvzoRz9ixowZjBw5kmXLlr237cEHH8zChQt56KGHGDVqFMOHD2fUqFE89dRTWxx37ty5DBs2jGHDhjF8+HDWrFlTsTY54ZuZNWPZsmVblEfu06cPtbW1jBs3jhtvvBGA5cuX89JLLzFixAj22Wcf5s2bx6JFi7jwwgs555xztjjuZZddxpVXXsnixYu599572W677SrSHnDCNzNrVkuVLyOCgw8+mJtuugmAG2+8kQkTJgCwevVqJkyYwODBgznjjDPe9y2g0f7778+ZZ57JtGnTWLVqVUVnyXLCNzNrxqBBg9i8ttcbb7zBiy++yMiRI+nbty9LlizhhhtueG+GqvPOO4/Ro0ezdOlSZs2axbp167Y47pQpU7j66qt5++232W+//XjyyScr0h5oY8KX9AFJQ4sKxsyss/j0pz/NW2+9xbXXXgskk5WcddZZnHLKKfTq1Yvjjz+eSy65hNWrVzNkyBAgucLffffdgWT6w+Y888wzDBkyhG9/+9vU19dXNOG3+l1C0hzgyHTbxcAKSXMj4syCYzMze0+WYZR5aiyP/I1vfIOpU6eyadOm95U7PvbYY5k0aRLnnXfee/t861vf4uSTT+byyy9nzJgxzR73iiuuYPbs2XTr1o2BAwdy2GGHVaQ9kKE8sqRFETFc0leAPSLi+5KWRETuV/ouj2zVxOWRS3N55C1Vojxyd0m7AV8AOvapATMzK1uWhH8hcA/wTEQskLQn8HSxYZmZWd5a7cOPiJuAm5q8fxb4fJFBmZlB65OCV5M8Zids9Qpf0l6S/iRpafp+qKTvtvvMZmYl1NTUsHLlylwSXVcXEaxcuZKampp2HSfLiP+rgLOBn6UnXiLpN8C/t+vMZmYl9O/fn4aGBlasWNHRoXQKNTU19O/fv13HyJLwe0XEQ5t9rdrQrrOambWiR48eDBgwoKPD2KpkuWn7qqSPAgEg6VhgeaFRmZlZ7rJc4Z8GTAf2kfS/wHPASYVGZWZmucsySudZ4DOStge2iYjMtTwlPQ+sATYCG7I+HGBmZvlrMeFLarZ0QmNffkRcnvEcoyPi1baHZmZmeSp1hd+7YlGYmVnhWkz4EXFBDscP4PeSAvhZREzffANJE4GJALW1tTmc0szMmpPlwas9Jc2StELSK5JuTcsrZLF/RHwCOAw4TdKBm28QEdMjoj4i6vv169fG8M3MLKsswzJ/A9wI7AZ8mKTMwnVZDh4RL6U/XwFuBj5ZXphmZtZeWRK+IuKXEbEhff2KdEx+yZ2k7SX1blwGxgJL2xeumZmVK8s4/NmSpgDXkyT644A7JO0MEBGvtbDfrsDN6aie7sBvIuLu9odsZmblyJLwj0t/fnWz9V8i+Q+g2f78dPz+vuWHZmZmecry4JWLWZiZbQWyzGnbDTgCqGu6fRsevDIzs04gS5fOLGAd8BiwqdhwzMysKFkSfv8iJiw3M7PKyjIs8y5JYwuPxMzMCpXlCv8BkuGV2wDrAQEREX0KjczMzHKVJeH/EPgU8Fh4ckkzsy4rS5fO08BSJ3szs64tyxX+cmCOpLuAdxpXelimmVnXkiXhP5e+eqYvMzPrgrI8aZtHXXwzM+tgWZ607Qd8CxgE1DSuj4gxBcZlZmY5y3LT9tfAk8AA4ALgeWBBgTGZmVkBsiT8vhHx38D6iJgbEV8C9is4LjMzy1mWm7br05/LJR0BvAT0Ly4kMzMrQpaE/++SdgTOAn4C9AHOKDQqMzPLXZZROreni6uB0QCSPDzTzKyLabUPX9IcSXVN3o/EN23NzLqcLF06FwF3S5oG7A4cDpxaaFRmZpa7LF0690j6GvAH4FVgeES8XHhkZmaWqyxdOueR3Kw9EDifpK7OEQXHZWZmOcvSpbML8MmIeBuYL+lu4GrgjkIjMzOzXGXp0pm02fsXfIVvZtb1tNilI+m+Jsu/3OzjBwqLyMzMClGqD3/7JsuDNvtMBcRiZmYFKpXwS81wlXn2K0ndJC2SdHvrW5uZWVFK9eHvJOlokv8UdpJ0TLpewI5tOMck4AmSkgxmZtZBSiX8ucCRTZbHN/lsXpaDS+oPHAH8B3BmOQGamVk+Wkz4EZHH07RXkEye0juHY5mZWTtkGYdfFknjgFciYqGkg0tsNxGYCFBbW1tUOFVj1rv3tb5RCeN7HpBTJGbW2WSZAKVc+wNHSnoeuB4YI+lXm28UEdMjoj4i6vv161dgOGZm1a2whB8R34mI/hFRBxwP/DkiTirqfGZmVlqmLh1Jo4C6pttHxLUFxWRmZgVoNeGnT9l+FFgMbExXB5A54UfEHGBO28MzM7O8ZLnCrwcGRkTmh63MzKzzydKHvxT4UNGBmJlZsbKWR35c0kPAO40rI+LIlncxM7POJkvCP7/oIMzMrHhZ6uHPlbQrMDJd9VBEvFJsWGZmlrcsUxx+AXgImAB8AXhQ0rFFB2ZmZvnK0qVzLjCy8apeUj/gj8DMIgMzM7N8ZRmls81mXTgrM+5nZmadSJYr/Lsl3QNcl74/DrizuJDMzKwIWW7ani3p8yTF0ARMj4ibC4/MzMxylamWTkT8FvhtwbGYmVmBWkz4ku6LiAMkreH9c9gKiIjwlIVmZl1IqRmvDkh/erYqM7OtQJZx+L/Mss7MzDq3LMMrBzV9I6k7MKKYcMzMrCgtJnxJ30n774dKeiN9rQH+DtxasQjNzCwXLSb8iLgo7b+/NCL6pK/eEdE3Ir5TwRjNzCwHWcbhf0fSB4CPATVN1s8rMjAzM8tXlikOvwJMAvqTTHO4HzAfGFNsaGZmlqcsN20nkZRGfiEiRgPDgRWFRmVmZrnLkvDXRcQ6AEnbRsSTwN7FhmVmZnnLUlqhQdJOwC3AHyS9DrxUbFhmZpa3LDdtj04Xz5c0G9gRuLvQqMzMLHelauns3Mzqx9KfOwCvFRKRmZkVotQV/kKSomkCaoHX0+WdgL8BAwqPzszMclPqwasBEbEncA8wPiJ2iYi+wDjgd5UK0MzM8pFllM7IiHhvhquIuAs4qLWdJNVIekjSo5KWSbqgPYGamVn7ZBml86qk7wK/IuniOYlkXtvWvAOMiYi1knoA90m6KyIeKD9cMzMrV5Yr/BOAfsDNJEMzP5iuKykSa9O3PdJXlNjFzMwKlGVY5mskT9u2maRuJDd//wm4MiIebGabicBEgNra2nJOY53MrHfva9f+43seUN5572/XaVs0flQxxzWrtFLDMq+IiMmSZtHMlXlEHNnawSNiIzAsfXDrZkmDI2LpZttMB6YD1NfX+xuAmVlBSl3hN85qdVl7TxIRqyTNAQ4FlrayuZmZFaDUnLYL059zyzmwpH7A+jTZbwd8BvjPsqI0M7N2K9Wl8xglbrJGxNBWjr0b8Iu0H38b4MaIuL2sKM3MrN1KdemMa8+BI2IJSSllMzPrBEp16bxQyUDMzKxYrY7Dl7SfpAWS1kp6V9JGSW9UIjgzM8tPlgevfkryoNXTwHbAV4CfFBmUmZnlL0tpBSLir5K6pePqfy6poEdczMysKFkS/luSegKLJV0CLAe2LzYsMzPLW5YunS+m2/0r8CawB/D5IoMyM7P8lRqH/2/ADU1G66wDXOLYzKyLKnWFvztwv6R5kr4uaZdKBWVmZvkrNePVGSRTG54HDAWWSLpL0j9L6l2pAM3MLB8l+/DTmvZzI+LrJH33VwBnAH+vRHBmZpafTMMyJQ0BjgeOI5nt6pwigzIzs/yVumn7MZIkfwKwEbgeGBsRz1YoNjMzy1GpK/x7gOuA4yLisQrFY2ZmBSlVPG3PSgZiZmbFyvLglZmZbQWc8M3MqkSLCV/Sn9KfnpbQzGwrUOqm7W6SDgKOlHQ9oKYfRsQjhUZmZma5KpXwvwdMAfoDl2/2WQBjigrKzMzyV2qUzkxgpqTzImJqBWMyM7MCtPqkbURMlXQkcGC6ak5E3F5sWGZmlrcsc9peBEwCHk9fk9J1ZmbWhWSppXMEMCwiNgFI+gWwCPhOkYGZmVm+so7D36nJ8o5FBGJmZsXKcoV/EbBI0mySoZkH4qt7M7MuJ8tN2+skzQFGkiT8b0fEy63tJ2kP4FrgQ8AmYHpE/Lh94ZqZWbky1cOPiOXAbW089gbgrIh4JJ0ha6GkP0TE420N0szM2q+wWjoRsbzxadyIWAM8QTJPrpmZdYCKFE+TVAcMBx6sxPnMzGxLJbt0JG0DLImIweWeQNIOwG+ByRHxRjOfTwQmAtTW1pZ7mk5h1rv3dXQI7da0DQueLOYcI/cp5riNFuyQ7+9hwZK27zNy7QG5xtAWs+7vsFO32fhRHR1B2xT1Z1upP4fWJjHfBDwqqaxMLKkHSbL/dUT8roVzTI+I+oio79evXzmnMTOzDLLctN0NWCbpIeDNxpURcWSpnSQJ+G/giYjYvPiamZlVWJaEf0GZx94f+CLwmKTF6bpzIuLOMo9nZmbtkGUc/lxJHwE+FhF/lNQL6JZhv/vYrIa+mZl1nCzF0/4FmAn8LF21O3BLkUGZmVn+sgzLPI2ke+YNgIh4GvhgkUGZmVn+siT8dyLi3cY3krqTzHhlZmZdSJaEP1fSOcB2kg4BbgJmFRuWmZnlLUvCnwKsAB4DvgrcCXy3yKDMzCx/WUbpbEonPXmQpCvnqYhwl46ZWRfTasKXdATwX8AzJMMsB0j6akTcVXRwZmaWnywPXv0QGB0RfwWQ9FHgDsAJ38ysC8nSh/9KY7JPPQu8UlA8ZmZWkBav8CUdky4uk3QncCNJH/4EYEEFYjMzsxyV6tIZ32T578BB6fIK4AOFRWRmZoVoMeFHxKmVDMTMzIqVZZTOAOB0oK7p9q2VRzYzs84lyyidW0jq2s8CNhUbjpmZFSVLwl8XEdMKj8TMzAqVJeH/WNL3gd8D7zSujIhHCovKzMxylyXhDyGZuWoM/+jSifS9mZl1EVkS/tHAnk1LJJuZWdeT5UnbR4Gdig7EzMyKleUKf1fgSUkLeH8fvodlmpl1IVkS/vcLj8LMzAqXpR7+3EoEYmZmxcrypO0a/jGHbU+gB/BmRPQpMjAzM8tXliv83k3fSzoK+GRhEZmZWSGyjNJ5n4i4BY/BNzPrcrJ06RzT5O02QD3/6OIxM7MuIssonaZ18TcAzwOfa20nSTOAcSQzZg0uKzozM8tNlj78cuviXwP8FLi2zP3NzCxHpaY4/F6J/SIippY6cETMk1RXZlxmZpazUlf4bzazbnvgy0BfoGTCz0rSRGAiQG1tbdnHmfXufXmEYxWw4MlWPqfM3+UO5e1mVi1KTXH4w8ZlSb2BScCpwPXAD1var60iYjowHaC+vt43g83MClKyD1/SzsCZwInAL4BPRMTrlQjMzMzyVaoP/1LgGJKr7yERsbZiUZmZWe5KPXh1FvBh4LvAS5LeSF9rJL3R2oElXQfMB/aW1CDpy/mEbGZm5SjVh9/mp3A32/+E9uxvZmb5aldSNzOzrsMJ38ysSjjhm5lVCSd8M7Mq4YRvZlYlnPDNzKqEE76ZWZVwwjczqxJO+GZmVcIJ38ysSjjhm5lVCSd8M7Mq4YRvZlYlnPDNzKqEE76ZWZVwwjczqxJO+GZmVcIJ38ysSjjhm5lVCSd8M7Mq4YRvZlYlnPDNzKqEE76ZWZVwwjczqxJO+GZmVaLQhC/pUElPSfqrpClFnsvMzEorLOFL6gZcCRwGDAROkDSwqPOZmVlpRV7hfxL4a0Q8GxHvAtcDnyvwfGZmVkKRCX934MUm7xvSdWZm1gG6F3hsNbMutthImghMTN+ulfRUgTFV0i7Aqx0dRIG25va5bV1TtbbtI1kPUmTCbwD2aPK+P/DS5htFxHRgeoFxdAhJD0dEfUfHUZStuX1uW9fktrWuyC6dBcDHJA2Q1BM4HritwPOZmVkJhV3hR8QGSf8K3AN0A2ZExLKizmdmZqUV2aVDRNwJ3FnkOTqxra6bajNbc/vctq7JbWuFIra4j2pmZlshl1YwM6sSTvjt1Fr5CElnSnpc0hJJf5KUeQhVR8taGkPSsZJCUpcZIZGlbZK+kP7ulkn6TaVjLFeGv5O1kmZLWpT+vTy8I+Ish6QZkl6RtLSFzyVpWtr2JZI+UekYy5WhbSembVoi6X5J+7b5JBHhV5kvkpvRzwB7Aj2BR4GBm20zGuiVLn8duKGj486rbel2vYF5wANAfUfHnePv7WPAIuAD6fsPdnTcObZtOvD1dHkg8HxHx92G9h0IfAJY2sLnhwN3kTwHtB/wYEfHnGPbRjX5+3hYOW3zFX77tFo+IiJmR8Rb6dsHSJ5H6AqylsaYClwCrKtkcO2UpW3/AlwZEa8DRMQrFY6xXFnaFkCfdHlHmnk+prOKiHnAayU2+RxwbSQeAHaStFtlomuf1toWEfc3/n2kzFzihN8+bS0f8WWSq4+uoNW2SRoO7BERt1cysBxk+b3tBewl6X8kPSDp0IpF1z5Z2nY+cJKkBpJRdKdXJrSKqJaSLmXlkkKHZVaBTOUjACSdBNQDBxUaUX5Ktk3SNsCPgFMqFVCOsvzeupN06xxMciV1r6TBEbGq4NjaK0vbTgCuiYgfSvoU8Mu0bZuKD69wmf9NdlWSRpMk/APauq+v8NsnU/kISZ8BzgWOjIh3KhRbe7XWtt7AYGCOpOdJ+ktv6yI3brP83hqAWyNifUQ8BzxF8h9AZ5elbV8GbgSIiPlADUmtlq1Bpn+TXZWkocDVwOciYmVb93fCb59Wy0ek3R4/I0n2XaUfGFppW0SsjohdIqIuIupI+hSPjIiHOybcNslS9uMWkhvuSNqFpIvn2YpGWZ4sbfsb8GkASR8nSfgrKhplcW4D/jkdrbMfsDoilnd0UHmQVAv8DvhiRPylnGO4S6cdooXyEZIuBB6OiNuAS4EdgJskAfwtIo7ssKAzyti2Lilj2+4Bxkp6HNgInF3OFVWlZWzbWcBVks4g6e44JdKhH52dpOtIutl2Se9BfB/oARAR/0VyT+Jw4K/AW8CpHRNp22Vo2/eAvsD/S3PJhmhjQTU/aWtmViXcpWNmViWc8M3MqoQTvplZlXDCNzOrEk74ZmZVwgnfOpSktQUcc6OkxZIelfSIpFF5nyNDDNdIei6NY7Gkb5ZxjGFdqZKldX4eh29bo7cjYhiApM8CF9ExJS3OjoiZ7dh/GEk5jsyzxknqHhEb2nFO24r5Ct86HUkfSecOaJxDoDZd/9G0kNkCSRdm/HbQB3g93X+H9HiPSHpM0ufS9VMlTWpy/v9ovCKXdHZ6viWSLkjXbS/pjvQbxFJJx7Whbf9f0sNpjf0LmqwfmdY4f1TSQ5J2BC4Ejku/IRwnaWdJt6SxPJA+Zo+k8yVNl/R74NqssVgV6uga0H5V9wtY28y6WcDJ6fKXgFvS5duBE9LlrzW3b/rZRmAx8CSwGhiRru8O9EmXdyF5GlNAHfBIun4bknryfYGxJLXjla6/naRm+eeBq5qcb8dmYrgGeC6NYzEwJF2/c/qzGzAHGEpSt/5ZYGT6WZ801lOAnzY55k+A76fLY4DF6fL5wEJgu47+ffrVuV/u0rHO6FPAMenyL0nq7TeuPypd/g1wWQv7N+3S+RRwraTBJIn7B5IOBDaRlM3dNSKel7QyrXu0K7AoIlZKGkuS9Belx92BpIDavcBlkv4TuD0i7m0hjua6dL4gaSJJQt+NZAKSAJZHxAKAiHgjjX3z4x1A8p8NEfFnSX3TbwIAt0XE2y3EYQa4D9+6hrLrf0TE/LT4WT+SGiv9SK7416dVPmvSTa8muaL+EDAjXSfgooj42ebHlTQiPd5Fkn4fERe2FoukAcC/kVzJvy7pmvT8ytjGUqV/38ywv1U59+FbZ3Q/SZVHgBOB+9LlB0ivcJt8XpKkfUi6T1aSzO70SprsRwNN5xe+GTgUGElSeIz055ck7ZAea3dJH5T0YeCtiPgVybeMrPOm9iFJzKsl7UoyTR0kXU8fljQyPU9vSd2BNSRlqBvNI/nzQNLBwKuN3wbMsvAVvnW0XmllwEaXA98EZkg6m6Rsb2PFw8nArySdBdxB0j/fnO0kLU6XRXI/YKOkXwOzJD3MP/r4AYiIdyXNBlZFxMZ03e/T8sHz0+6VtcBJwD8Bl0raBKwnmau4VRHxqKRFwDKSPvv/aXLu44CfSNoOeBv4DDAbmJK25SKSvvqfS1pCUgny5CznNWvkapnWZUjqRdI/H5KOJ7mB29w8u+UcexvgEWBCRDydxzHNOhtf4VtXMgL4qZLL7VUkI3jaTdJAkhE4NzvZ29bMV/hmZlXCN23NzKqEE76ZWZVwwjczqxJO+GZmVcIJ38ysSjjhm5lVif8D9pWbJS8KGz4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(lbf_circs, alpha = 0.3)\n",
    "plt.hist(lbf_ovals, alpha = 0.3)\n",
    "plt.legend(['Circles','Ovals'])\n",
    "plt.xlabel('Log Bayes Factor')\n",
    "plt.ylabel('Number of Validation Examples')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On average, we notice higher Bayes factors for the circles, which was to be expected. To quantitavely assess the Bayes factor classifier, we can compute the AUC. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9125"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import metrics\n",
    "\n",
    "class_labels = [0]*data_sz + [1]*data_sz\n",
    "pred = lbf_ovals + lbf_circs\n",
    "fpr, tpr, thresholds = metrics.roc_curve(class_labels, pred, pos_label=1)\n",
    "metrics.auc(fpr,tpr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The AUC is close to 1, which indicates good performance of our classifier. "
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}