08-Mapper.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapper\n",
    "\n",
    "In this notebook, we will explore the Mapper algorithm for graph-based dimension reduction and data visualization. A nice implementation of the Mapper algorithm is built into the `giotto-tda` library. The first part of the notebook is adapted from a tutorial on the `giotto` [Github](https://github.com/giotto-ai/giotto-tda). As was mentioned last week, they have several other nice notebooks describing applications of TDA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from gtda.plotting import plot_point_cloud\n",
    "\n",
    "from gtda.mapper import (\n",
    "    CubicalCover,\n",
    "    OneDimensionalCover,\n",
    "    FirstSimpleGap,\n",
    "    Eccentricity,\n",
    "    Entropy,\n",
    "    FirstHistogramGap,\n",
    "    make_mapper_pipeline,\n",
    "    Projection,\n",
    "    plot_static_mapper_graph,\n",
    "    plot_interactive_mapper_graph\n",
    ")\n",
    "\n",
    "from sklearn import datasets\n",
    "from sklearn.cluster import DBSCAN\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Toy Data\n",
    "\n",
    "We'll start with a basic example which will allow us to easily interpret the results of the mapper algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data, _ = datasets.make_circles(n_samples=5000, noise=0.05, factor=0.3, random_state=42)\n",
    "\n",
    "plot_point_cloud(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Mapper Algorithm\n",
    "\n",
    "Recall that a Mapper visualization of a finite metric space $(X,d_X)$ (such as a point cloud in $\\mathbb{R}^n$) requires several parameters:\n",
    "\n",
    "- A reference topological space $Y$ --- typically $\\mathbb{R}$\n",
    "\n",
    "- A function $f:X \\to Y$ --- e.g., a data-driven function, a density estimate, eccentricity, etc.\n",
    "\n",
    "- An open cover of $Y$ --- we'll generally pick from a built-in family of covers of $\\mathbb{R}$\n",
    "\n",
    "- A clustering algorithm --- our definition in class always used clusters coming from a Vietoris-Rips complex at a chosen scale $\\epsilon > 0$; as was mentioned in class, we can choose from other clustering algorithms as well.\n",
    "\n",
    "For each set $V_i$ in the open cover of $Y$, we pull back to an open set $U_i = f^{-1}(V_i)$ in $Y$. The open set is divided into clusters $U_i^{(1)},\\ldots,U_i^{(n_i)}$ via the clustering algorithm. The Mapper graph is the 1-skeleton of the *nerve* of the collection $\\{U_i^{(j)}\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's choose various parameters now."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter Function\n",
    "\n",
    "We can choos any `scikit-learn` [transformer](https://scikit-learn.org/stable/data_transforms.html). There are also some functions built into `giotto-tda`---see the [documentation](https://giotto-ai.github.io/gtda-docs/latest/modules/mapper.html). \n",
    "\n",
    "Here we are using a function $f:X \\to \\mathbb{R}$ which projects onto the second coordinate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_func = Projection(columns=[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cover\n",
    "\n",
    "We'll choose from a built-in family of covers. Here we are just choosing a number of intervals and how much they overlap (similar to the picture we drew in lecture)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cover = OneDimensionalCover(n_intervals=10, overlap_frac=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Clustering\n",
    "\n",
    "The default clustering algorithm is called `DBSCAN`---this is more refined than the Vietoris-Rips-based clustering described in class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusterer = DBSCAN()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mapper Pipeline\n",
    "\n",
    "Now we create a Mapper pipeline which takes in all of these parameters. There is also an `n_jobs` parameter, in case we wanted to run some things in parallel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipe = make_mapper_pipeline(\n",
    "    filter_func=filter_func,\n",
    "    cover=cover,\n",
    "    clusterer=clusterer,\n",
    "    verbose=False,\n",
    "    n_jobs=1,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualise the Mapper graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a plot of the resulting mapper graph with default options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot_static_mapper_graph(pipe, data)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these choices of parameters aren't capturing the correct topology. The Mapper graph *is* observing the two holes in the data, but it's not seeing the two connected components. Let's play around with parameters to see if we can get a Mapper graph that looks closer to the original data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we can change around the number of intervals and overlap fraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_func = Projection(columns=[1])\n",
    "cover = OneDimensionalCover(n_intervals=20, overlap_frac=0.1)\n",
    "clusterer = DBSCAN()\n",
    "\n",
    "pipe = make_mapper_pipeline(\n",
    "    filter_func=filter_func,\n",
    "    cover=cover,\n",
    "    clusterer=clusterer,\n",
    "    verbose=False,\n",
    "    n_jobs=1,\n",
    ")\n",
    "\n",
    "fig = plot_static_mapper_graph(pipe, data)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I can't find parameters that will get the above to work. \n",
    "\n",
    "We can also change the clustering algorithm. Here we use the `FirstSimpleGap` function. This uses the Vietoris-Rips construction, but automatically chooses $\\epsilon$ to be the height of the first 'significant gap' in the hierarchical clustering dendrogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "filter_func = Projection(columns=[1])\n",
    "cover = OneDimensionalCover(n_intervals=30, overlap_frac=0.3)\n",
    "clusterer = FirstSimpleGap()\n",
    "\n",
    "pipe = make_mapper_pipeline(\n",
    "    filter_func=filter_func,\n",
    "    cover=cover,\n",
    "    clusterer=clusterer,\n",
    "    verbose=False,\n",
    "    n_jobs=1,\n",
    ")\n",
    "\n",
    "fig = plot_static_mapper_graph(pipe, data)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This seems to give a correct result!\n",
    "\n",
    "We can also change the filter function to, say, $p$-eccentricity (as introduced in class). This picks up connected components, but not cycles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_func = Eccentricity(exponent = 2)\n",
    "cover = OneDimensionalCover(n_intervals=30, overlap_frac=0.1)\n",
    "clusterer = FirstSimpleGap()\n",
    "\n",
    "pipe = make_mapper_pipeline(\n",
    "    filter_func=filter_func,\n",
    "    cover=cover,\n",
    "    clusterer=clusterer,\n",
    "    verbose=False,\n",
    "    n_jobs=1,\n",
    ")\n",
    "\n",
    "fig = plot_static_mapper_graph(pipe, data)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we could change the target space $Y$ *and* the function $f:X \\to Y$. We'll use $Y = \\mathbb{R}^2$ and $f$ the identity function---this is not realistic for high-dimensional data, but we can use it for this low-dimensional example.\n",
    "\n",
    "Since we are changing the target space, we also need to change the type of cover. We'll cover $Y$ by squares with a prescribed overlap percentage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_func = Projection(columns=[0,1])\n",
    "cover = CubicalCover(n_intervals=10, overlap_frac=0.2)\n",
    "clusterer = DBSCAN()\n",
    "\n",
    "pipe = make_mapper_pipeline(\n",
    "    filter_func=filter_func,\n",
    "    cover=cover,\n",
    "    clusterer=clusterer,\n",
    "    verbose=False,\n",
    "    n_jobs=1,\n",
    ")\n",
    "\n",
    "fig = plot_static_mapper_graph(pipe, data)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This produces a very good result, which makes sense since we are using a filtration which faithfully represents the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization Parameters\n",
    "\n",
    "There are lots of options built into `giotto-tda` for visualizing Mapper graphs. We'll explore some options here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Node Coloring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotly_params = {\"node_trace\": {\"marker_colorscale\": \"Blues\"}}\n",
    "fig = plot_static_mapper_graph(\n",
    "    pipe, data, color_by_columns_dropdown=True, plotly_params=plotly_params\n",
    ")\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Graph Layout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot_static_mapper_graph(\n",
    "    pipe, data, layout=\"fruchterman_reingold\", color_by_columns_dropdown=True\n",
    ")\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting Dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot_static_mapper_graph(pipe, data, layout_dim=3, color_by_columns_dropdown=True)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Node scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_scale = 50\n",
    "fig = plot_static_mapper_graph(pipe, data, layout_dim=3, node_scale=node_scale)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the Mapper pipeline\n",
    "\n",
    "Behind the scenes of ``plot_static_mapper_graph`` is a ``MapperPipeline`` object ``pipe`` that can be used like a typical ``scikit-learn`` estimator. For example, to extract the underlying graph data structure we can do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph = pipe.fit_transform(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting graph is a [python-igraph](https://igraph.org/python/) object which stores node metadata in the form of attributes. We can access this data as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graph.vs.attributes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here ``'pullback_set_label'`` and ``'partial_cluster_label'`` refer to the interval and cluster sets described above. ``'node_elements'`` refers to the indices of our original data that belong to each node. For example, to find which points belong to the first node of the graph we can access the desired data as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_id = 0\n",
    "node_elements = graph.vs[\"node_elements\"]\n",
    "\n",
    "print(f\"\"\"\n",
    "Node ID: {node_id}\n",
    "Node elements: {node_elements[node_id]}\n",
    "Data points: {data[node_elements[node_id]]}\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating custom filter functions\n",
    "\n",
    "We can create custom filter functions. These will act row-wise on the input data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_func = np.sum\n",
    "\n",
    "pipe = make_mapper_pipeline(\n",
    "    filter_func=filter_func,\n",
    "    cover=cover,\n",
    "    clusterer=clusterer,\n",
    "    verbose=False,\n",
    "    n_jobs=1,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plot_static_mapper_graph(pipe, data)\n",
    "fig.show(config={'scrollZoom': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Visualization\n",
    "\n",
    "A useful feature is live visualization, which can be used to tune parameters and see the result in real time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipe = make_mapper_pipeline()\n",
    "\n",
    "# Generate interactive widget\n",
    "plot_interactive_mapper_graph(pipe, data, color_by_columns_dropdown=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing Mapper on Some Datasets\n",
    "\n",
    "Next let's try to produce Mapper graphs for some real data. The `sklearn` package has several classic toy datasets that can be loaded. The available datasets can be found [here](https://scikit-learn.org/stable/datasets/toy_dataset.html).\n",
    "\n",
    "We'll start with the `iris` dataset which contains measurements of flower biology coming from three classes of iris'.\n",
    "\n",
    "The syntax for loading the data looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "iris = load_iris()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This loads a dictionary, and we want to collect the data matrix and labels for analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = load_iris()\n",
    "iris"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = iris['data']\n",
    "y = iris['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the shape of the data matrix, we see that there are 150 samples, with 4 features per sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 3 classes, indexed by the labels in `y`. The distribution of classes is even."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(y)\n",
    "plt.hist(y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the Mapper graph and play with parameters. Note that we can color the nodes by ground truth label, rather than function value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filter_func = Projection(columns = [0])\n",
    "clusterer = FirstSimpleGap()\n",
    "\n",
    "pipe = make_mapper_pipeline(filter_func = filter_func, clusterer = clusterer)\n",
    "plot_interactive_mapper_graph(pipe, X, color_by_columns_dropdown=True, color_variable = y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot 2D projections of the data. We see that (for certain choices of parameters), the mapper graph does seem to reflect the shape of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj_coord1 = 0\n",
    "proj_coord2 = 1\n",
    "\n",
    "plt.scatter(X[:,proj_coord1],X[:,proj_coord2],c = y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Group Work: Other Datasets\n",
    "\n",
    "In breakout groups, try exploring some other datasets. You can either load one of the other `sklearn` datasets, or use a dataset of your own that you'd like to test this on.\n",
    "\n",
    "Feel free to work independently in your groups and compare results, or all work together on a single dataset. I will bounce between rooms and try to help troubleshoot as necessary.\n",
    "\n",
    "Probably the most familiar of the datasets is `digits`, containing low res samples from the famous MNIST handwritten digit dataset.\n",
    "\n",
    "**Hint:** It can be useful to do some simple linear dimension reduction via PCA before applying the topological dimension reduction. This step is not necessary, but I had some luck with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra Time: Exploring Neural Network Activations\n",
    "\n",
    "There is a cool application of Mapper in [this paper](http://www.sci.utah.edu/~beiwang/publications/TopoAct_BeiWang_2021.pdf) by Rathore et al. They use Mapper graphs to visualize the structural changes that occur as data is passed through the layers of a trained neural network. The paper comes with an interactive demo, available [here](https://tdavislab.github.io/TopoAct/). If you have extra time, explore this webtool with your group---try to understand what the Mapper graphs represent and see if you can find any interesting features as you explore the neural network."
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}