Approximate Sparse Filtrations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximate Sparse Filtrations\n",
    "\n",
    "In this module, we will explore an approximation algorithm which is meant to reduce the run time of the persistence algorithm [1].  For more information on this algorithm, please view the following video:\n",
    "\n",
    "https://www.youtube.com/watch?v=3WxuSwQhAgU\n",
    "\n",
    "\n",
    "[1] Cavanna, Nicholas J., Mahmoodreza Jahanseir, and Donald R. Sheehy. \"<a href = \"https://arxiv.org/pdf/1506.03797.pdf\">A geometric perspective on sparse filtrations</a>.\" Proceedings of the Canadian Conference on Computational Geometry (CCCG 2015)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ripser import ripser\n",
    "from persim import plot_diagrams, wasserstein, wasserstein_matching\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "from scipy import sparse\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define a \"greedy permutation,\" or a function which performs furthest points sampling, a key step in the algorithm used to determine \"insertion radii\" $\\lambda_i$ for each point (this is similar to a cover tree).  For an animation of how this works, please visit: \n",
    "\n",
    "https://gist.github.com/ctralie/128cc07da67f1d2e10ea470ee2d23fe8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getGreedyPerm(D):\n",
    "    \"\"\"\n",
    "    A Naive O(N^2) algorithm to do furthest points sampling\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    D : ndarray (N, N) \n",
    "        An NxN distance matrix for points\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    lamdas: list\n",
    "        Insertion radii of all points\n",
    "    \"\"\"\n",
    "    \n",
    "    N = D.shape[0]\n",
    "    #By default, takes the first point in the permutation to be the\n",
    "    #first point in the point cloud, but could be random\n",
    "    perm = np.zeros(N, dtype=np.int64)\n",
    "    lambdas = np.zeros(N)\n",
    "    ds = D[0, :]\n",
    "    for i in range(1, N):\n",
    "        idx = np.argmax(ds)\n",
    "        perm[i] = idx\n",
    "        lambdas[i] = ds[idx]\n",
    "        ds = np.minimum(ds, D[idx, :])\n",
    "    return lambdas[perm]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define a function which, given a distance matrix representing a point cloud, a set of insertion radii, and an approximation factor $\\epsilon$, returns a sparse distance matrix with re-weighted edges, whose persistence diagrams are each guaranteed to be a $(1+\\epsilon)$ multiplicative approximation of the true persistence diagrams (see [1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getApproxSparseDM(lambdas, eps, D):\n",
    "    \"\"\"\n",
    "    Purpose: To return the sparse edge list with the warped distances, sorted by weight\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lambdas: list\n",
    "        insertion radii for points\n",
    "    eps: float\n",
    "        epsilon approximation constant\n",
    "    D: ndarray\n",
    "        NxN distance matrix, okay to modify because last time it's used\n",
    "    \n",
    "    Return\n",
    "    ------\n",
    "    DSparse: scipy.sparse\n",
    "        A sparse NxN matrix with the reweighted edges\n",
    "    \"\"\"\n",
    "    N = D.shape[0]\n",
    "    E0 = (1+eps)/eps\n",
    "    E1 = (1+eps)**2/eps\n",
    "    \n",
    "    #Create initial sparse list candidates (Lemma 6)\n",
    "    nBounds = ((eps**2+3*eps+2)/eps)*lambdas #Search neighborhoods    \n",
    "    D[D > nBounds[:, None]] = np.inf #Set all distances outside of search neighborhood to infinity\n",
    "    [I, J] = np.meshgrid(np.arange(N), np.arange(N))\n",
    "    idx = I < J\n",
    "    I = I[(D < np.inf)*(idx == 1)]\n",
    "    J = J[(D < np.inf)*(idx == 1)]\n",
    "    D = D[(D < np.inf)*(idx == 1)]\n",
    "    \n",
    "    #Prune sparse list and update warped edge lengths (Algorithm 3 pg. 14)\n",
    "    minlam = np.minimum(lambdas[I], lambdas[J])\n",
    "    maxlam = np.maximum(lambdas[I], lambdas[J])\n",
    "    #Rule out edges between vertices whose balls stop growing before they touch\n",
    "    #or where one of them would have been deleted.  M stores which of these\n",
    "    #happens first\n",
    "    M = np.minimum((E0 + E1)*minlam, E0*(minlam + maxlam))\n",
    "    #M = E0*(minlam+maxlam)\n",
    "    \n",
    "    t = np.arange(len(I))\n",
    "    t = t[D <= M]\n",
    "    (I, J, D) = (I[t], J[t], D[t])\n",
    "    minlam = minlam[t]\n",
    "    maxlam = maxlam[t]\n",
    "    \n",
    "    #Now figure out the metric of the edges that are actually added\n",
    "    t = np.ones(len(I))\n",
    "    t[D <= 2*minlam*E0] = 0 #If cones haven't turned into cylinders, metric is unchanged\n",
    "    #Otherwise, if they meet before the M condition above, the metric is warped\n",
    "    D[t == 1] = 2.0*(D[t == 1] - minlam[t == 1]*E0) #Multiply by 2 convention\n",
    "    return sparse.coo_matrix((D, (I, J)), shape=(N, N)).tocsr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's set up a point cloud we can test this on, which has enough points for ripser to start slowing down a bit.  We'll perform the full rips filtration on this point cloud as a ground truth, and we will time it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NPoints = 2000\n",
    "t = np.linspace(0, 2*np.pi, NPoints+1)[0:NPoints]\n",
    "X = np.zeros((NPoints, 2))\n",
    "X[:, 0] = np.cos(t)\n",
    "X[:, 1] = np.sin(2*t)\n",
    "X += 0.1*np.random.randn(NPoints, 2)\n",
    "plt.scatter(X[:, 0], X[:, 1])\n",
    "\n",
    "\n",
    "tic = time.time()\n",
    "resultfull = ripser(X)\n",
    "toc = time.time()\n",
    "timefull = toc-tic\n",
    "print(\"Elapsed Time: %.3g seconds, %i Edges added\"%(timefull, resultfull['num_edges']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's run an approximate version and plot H1 for both next to each other\n",
    "\n",
    "#### Questions\n",
    "* What happens to the approximation when you change $\\epsilon$?  What happens to the runtime?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps = 0.05\n",
    "\n",
    "\n",
    "#Cmpute the sparse filtration\n",
    "tic = time.time()\n",
    "#First compute all pairwise distances and do furthest point sampling\n",
    "D = pairwise_distances(X, metric='euclidean')\n",
    "lambdas = getGreedyPerm(D)\n",
    "#Now compute the sparse distance matrix\n",
    "DSparse = getApproxSparseDM(lambdas, eps, D)\n",
    "#Finally, compute the filtration\n",
    "resultsparse = ripser(DSparse, distance_matrix=True)\n",
    "toc = time.time()\n",
    "timesparse = toc-tic\n",
    "percent_added = 100.0*float(resultsparse['num_edges'])/resultfull['num_edges']\n",
    "print(\"Elapsed Time: %.3g seconds, %i Edges added\"%(timesparse, resultsparse['num_edges']))\n",
    "\n",
    "\n",
    "#Plot the sparse distance matrix and edges that were added\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.subplot(121)\n",
    "D = pairwise_distances(X, metric='euclidean')\n",
    "plt.imshow(D)\n",
    "plt.title(\"Original Distance Matrix: %i Edges\"%resultfull['num_edges'])\n",
    "plt.subplot(122)\n",
    "DSparse = DSparse.toarray()\n",
    "DSparse = DSparse + DSparse.T\n",
    "plt.imshow(DSparse)\n",
    "plt.title(\"Sparse Distance Matrix: %i Edges\"%resultsparse['num_edges'])\n",
    "\n",
    "\n",
    "\n",
    "#And plot the persistence diagrams on top of each other\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.subplot(121)\n",
    "plot_diagrams(resultfull['dgms'], show=False)\n",
    "plt.title(\"Full Filtration: Elapsed Time %g Seconds\"%timefull)\n",
    "plt.subplot(122)\n",
    "plt.title(\"Sparse Filtration (%.3g%% Added)\\nElapsed Time %g Seconds\"%(percent_added, timesparse))\n",
    "plot_diagrams(resultsparse['dgms'], show=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compute the Wasserstein distance between the two diagrams to see how close they are.  Try this with different approximations above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Show the Wasserstein distance for H1\n",
    "I1 = resultfull['dgms'][1]\n",
    "I2 = resultsparse['dgms'][1]\n",
    "matchdist, (matchidx, D) = wasserstein(I1, I2, matching=True)\n",
    "plt.figure(figsize=(5, 5))\n",
    "#Plot wasserstein matching\n",
    "wasserstein_matching(I1, I2, matchidx, labels = ['exact', '$\\epsilon = %.3g$'%eps])\n",
    "plt.title(\"Wasserstein Dist: %.3g\"%matchdist)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}