day-4.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 4: Applications of persistent homology\n",
    "\n",
    "Today, we will explore some simple applications of persistent homology"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing distance between persistence diagrams\n",
    "\n",
    "First, lets compute the bottleneck and Wasserstein distance between persistence diagrams. Use Gudhi to write two functions - one for computing bottleneck distance and another for computing Wasserstein distance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy\n",
    "import gudhi\n",
    "import ot \n",
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_hdf5_data(filename):\n",
    "    dataset = None\n",
    "    labels = None\n",
    "    with h5py.File(filename,'r') as h5f: \n",
    "        dataset = h5f['data'][:]\n",
    "        labels = h5f['label'][:]\n",
    "    dataset = dataset - np.expand_dims(np.mean(dataset, axis=0), 0)  # center\n",
    "    return dataset, labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_bottleneck(dgm1, dgm2):\n",
    "    ''' Computes the bottleneck distance between two persistence diagrams\n",
    "\n",
    "    Parameters:\n",
    "        dgm1, dgm2: (N, 2) numpy array representing persistence points\n",
    "\n",
    "    Returns:\n",
    "        bottleneck distance \n",
    "\n",
    "    '''\n",
    "\n",
    "    raise NotImplementedError('Please implement')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_wasserstein(dgm1, dgm2):\n",
    "    ''' Computes the Wasserstein distance between two persistence diagrams\n",
    "\n",
    "    Parameters:\n",
    "        dgm1, dgm2: (N, 2) numpy array representing persistence points\n",
    "\n",
    "    Returns:\n",
    "        Wasserstein distance \n",
    "\n",
    "    '''\n",
    "    raise NotImplementedError('Please implement')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: try computing bottleneck and Wasserstein distance between a couple of persistence diagrams. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\pi: \\mathbb{R}^2 \\to \\mathbb{R}^2$ be the projection of any point onto the diagonal. Then, given two persistence diagrams, $P$ and $Q$:\n",
    "\n",
    "(1) show that the Wasserstein distance between $P \\cup \\pi(Q)$ and $Q \\cup \\pi(P)$ is the same the Wasserstein distance between $P$ and $Q$. \n",
    "\n",
    "(2) Use [Python OT](https://pythonot.github.io/index.html) to implement this alternate Wasserstein distance computation between $P$ and $Q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def projection_to_diagonal(P):\n",
    "    ''' Computes projection of a point set in \\R^2 to the diagonal \n",
    "\n",
    "    Parameters:\n",
    "        P: point set in \\R^2 formatted as (N, 2) np.array \n",
    "\n",
    "    Returns:\n",
    "        (N, 2) np.array of points on diagonal\n",
    "\n",
    "    '''\n",
    "    raise NotImplementedError('Please implement')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_wasserstein_alternate(dgm1, dgm2):\n",
    "    ''' Computes the Wasserstein distance between two persistence diagrams\n",
    "\n",
    "    Parameters:\n",
    "        dgm1, dgm2: (N, 2) numpy array representing persistence points\n",
    "\n",
    "    Returns:\n",
    "        Wasserstein distance \n",
    "\n",
    "    '''\n",
    "    raise NotImplementedError('Please implement')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point cloud classification using persistence diagrams \n",
    "\n",
    "Now, let's try doing some point cloud classification using persistence diagrams. We can do this in three steps: \n",
    "(1) we will compute persistence diagrams for our point cloud dataset. \n",
    "(2) we will compute the bottleneck/Wasserstein distance between each pair of persistence diagrams.   \n",
    "(3) we will use the sklearn k-nearest neighbor classifier to do point cloud classification. \n",
    "    \n",
    "We will be working with the ModelNet40 point cloud dataset. The training point clouds/labelsare in the `day3/train` folder and the testing point clouds/labels are in the `day3/test` folder. Use the `load_hdf5` function to load the point cloud data. Additionally, I've written the classifier already for you so all you need to do is step 1 and 2 and input the computed distance matrices in the correct places. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "k = 3\n",
    "# TODO: compute X_train, should be of shape (N_train, N_train)\n",
    "X_train = None\n",
    "# TODO: compute Y_train, should be shape (N_train, )\n",
    "Y_train = None\n",
    "\n",
    "# TODO: compute X_test, should be of shape (N_test, N_train)\n",
    "X_test = None\n",
    "# TODO: compute Y_test, should be shape (N_test, )\n",
    "Y_test = None\n",
    "\n",
    "neigh = KNeighborsClassifier(n_neighbors=k, metric='precomputed')\n",
    "\n",
    "neigh.fit(X_train, Y_train)\n",
    "\n",
    "predictions = neigh.predict(X_test)\n",
    "\n",
    "print(\"Classification accuracy:\", accuracy_score(Y_test, predictions))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have finished with a k-nearest neighbor classifier, try doing classification using an SVM from sklearn. If you've finished both, try using both methods with a image data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import svm\n",
    "\n",
    "k = 3\n",
    "# TODO: compute X_train, should be of shape (N_train, N_train)\n",
    "X_train = None\n",
    "# TODO: compute Y_train, should be shape (N_train, )\n",
    "Y_train = None\n",
    "\n",
    "# TODO: compute X_test, should be of shape (N_test, N_train)\n",
    "X_test = None\n",
    "# TODO: compute Y_test, should be shape (N_test, )\n",
    "Y_test = None\n",
    "\n",
    "clf = svm.SVC()\n",
    "\n",
    "clf.fit(X_train, Y_train)\n",
    "\n",
    "predictions = clf.predict(X_test)\n",
    "\n",
    "print(\"Classification accuracy:\", accuracy_score(Y_test, predictions))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "msri",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "d53e343ab014b81a4b593937664e6856782a01e1493065db0a4fc92ef9a58926"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}