exercises/Week-3-density-based-clustering-23.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 3 - Density Based Clustering\n",
    "\n",
    "During this weeks exercises, we will be working with density based clustering.\n",
    "In particular, we will be working with the algorithms DBSCAN, DENCLUE.\n",
    "\n",
    "First, lets import some stuff and plot the data that we are going to use today."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Local imports\n",
    "import sys\n",
    "sys.path.append('.')\n",
    "from utilities.load_data import load_iris_PC, load_t7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Theoretical questions\n",
    "- Please provide a brief description of what characterises density-based clustering as a clustering approach? How do both DBSCAN and DENCLUE define clusters, and what is the core difference between the two? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2: DBSCAN in Action\n",
    "Consider Figure 15.12 and answer the following questions, assuming that we use the\n",
    "Euclidean distance between points, and that $\\epsilon=2$ and $minpts = 3$\n",
    "1. List all the core points.\n",
    "1. Is $a$ directly density reachable from $d$?\n",
    "1. Is $o$ density reachable from $i$? Show the intermediate points on the chain or the point where the chain breaks.\n",
    "1. Is density reachable a symmetric relationship, that is, if $x$ is density reachable from $y$, does it imply that $y$ is density reachable from $x$? Why or why not?\n",
    "1. Is $l$ density connected to $x$? Show the intermediate points that make them density connected or violate the property, respectively.\n",
    "1. Is density connected a symmetric relationship?\n",
    "1. Show the density-based clusters and the noise points.\n",
    "1. If we use the manhattan distance instead, what is then the core points?\n",
    "\n",
    "<img src=\"graphics/15.12.png\" width=\"501\" />\n",
    "\n",
    "We have included the points in the code below, if you want to use `numpy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array([\n",
    "    [5., 10., 11., 6., 10., 12., 13., 5., 10., 13., 6., 9., 11., 14., 15., 2., 3., 5., 6., 7., 15., 3., 7., 8.],\n",
    "    [8.,  8.,  8., 7.,  7.,  7.,  7., 6.,  6.,  6., 5., 4.,  5.,  6.,  5., 4., 4., 4., 4., 4.,  4., 3., 3., 2.]\n",
    "]).T\n",
    "labels = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3: Implement DBSCAN\n",
    "\n",
    "In this exercise, we will try to implement the DBSCAN algorithm. \n",
    "You are allowed to structure your code however you want. \n",
    "The code below is inspired by [Zaki, p. 377] and serves as inspiration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def L2_norm(x, y):\n",
    "    return np.sqrt(np.sum((x-y)**2))\n",
    "\n",
    "def densityConnected(x, k, ...): # You might need more parameters here.\n",
    "    for neighbor in x.Nx:\n",
    "        # Connect (potentially recursively)\n",
    "        pass\n",
    "\n",
    "def dbscan(X, e, m, dist_fn=L2_norm):\n",
    "    \"\"\"\n",
    "        Args:\n",
    "            X:       Matrix of shape [n, d] with data points on the rows.\n",
    "            e:       Epsilon distance for neighborhood calculations.\n",
    "            m:       Minimum neighbors in epsilon neighborhood for a point to be a core point.\n",
    "            dist_fn: Distance function to be used.\n",
    "        \n",
    "        Returns:\n",
    "            clusters:   A vector of shape [n,] with integers, indicating cluster assignments.\n",
    "                        Let clusters[i] == -1 if point x_i is an outlier and a non-negative \n",
    "                        integer corresponding to the cluster index of point x_i otherwise.\n",
    "    \"\"\"\n",
    "    \n",
    "    # TODO, code here.\n",
    "    \n",
    "    core = []\n",
    "    for x in X:\n",
    "        # Compute neighborhoods and identify cores\n",
    "        pass\n",
    "\n",
    "    k = 0\n",
    "    for c in core:\n",
    "        # Assign unassigned cores to clusters\n",
    "        # using densityConnect\n",
    "        pass\n",
    "\n",
    "    clusters = []\n",
    "    for i in range(k):\n",
    "        # Build cluster indicator vector\n",
    "        pass\n",
    "    \n",
    "    # End TODO\n",
    "    \n",
    "    return clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is some code to try out your implementation.\n",
    "You can, of course, fiddle with the parameters and see what happens to the clusters.\n",
    "The parameters provided below whould work relatively well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test clustering on iris_PC dataset\n",
    "# clustering = dbscan(X, 0.5, 5) # Works for iris_2PC dataset\n",
    "slow = False\n",
    "if slow:\n",
    "    X, _ = load_t7()\n",
    "    clustering = dbscan(X, e=15, m=10)\n",
    "else: \n",
    "    X, _ = load_iris_PC()\n",
    "    clustering = dbscan(X, e=0.4, m=5)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.scatter(*(X.T), c=clustering)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4: DENCLUE Calculations\n",
    "This exercise relates to the DENCLUE algorithm introduced in Section 15.3 in [Zaki].\n",
    "Consider the points shown in Figure 15.13. For the Gaussian kernel\n",
    "\n",
    "$$\n",
    "K(\\mathbf{z})=\\frac{1}{(2 \\pi)^{d / 2}} \\exp \\left\\{-\\frac{\\mathbf{z}^{T} \\mathbf{z}}{2}\\right\\},\n",
    "$$\n",
    "\n",
    "<img src=\"graphics/15.13.png\" width=\"400\" />\n",
    "\n",
    "answer the following questions assuming that $h = 2$:\n",
    "\n",
    "1. What is the probability density at e?\n",
    "2. What is the gradient at e? (Try to actually derive the gradient of $\\hat f(x)$ your self)\n",
    "3. List all the density attractors for this dataset.\n",
    "\n",
    "As before, if you want to use `numpy`, we have included the points below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array([\n",
    "    [5., 6., 6., 2., 3., 5., 7., 9., 3., 8., 7.],\n",
    "    [8., 7., 5., 4., 4., 4., 4., 4., 3., 2., 5.]\n",
    "]).T\n",
    "labels = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5: Contingency table \n",
    "\n",
    "Provide an implementation of the contingency table and plot the table varying the DBScan parameters. Use iris_pc dataset. \n",
    "1. Can you use the contingengy table to evaluate the quality of clustering? \n",
    "2. Can you compute purity based on the contingency table?\n",
    "3. Can you find the best parameter setting using purity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X,ground_truth = load_iris_PC()\n",
    "n=ground_truth.shape[0]\n",
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "def convert(Y): # convert nxc one-hot matrix to nx1 vector\n",
    "    res=np.ndarray(Y.shape[0])\n",
    "    for i in range(Y.shape[0]):\n",
    "        res[i]=np.where(Y[i]==1)[0]\n",
    "    return res\n",
    "\n",
    "ground_truth=convert(ground_truth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contingency_table(C, T): \n",
    "    \"\"\"\n",
    "    Args:\n",
    "        C:       Clusters obtained by a clustering algorithm as a nx1 vector\n",
    "        T:       Ground-truth cluster assignments as a nx1 vector\n",
    "\n",
    "    Returns:\n",
    "        ctable:   a num_clusters_of_C x num_cluster_of_T matrix containing the overlaps among the different clusters\n",
    "    \"\"\"\n",
    "    \n",
    "    ctable = None\n",
    "    #YOUR CODE HERE\n",
    "    #YOUR CODE HERE\n",
    "    return ctable\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def purity(ctable,n):\n",
    "    \"\"\"\n",
    "    Args:\n",
    "        ctable:   a num_clusters_of_C x num_cluster_of_T matrix containing the overlaps among the different clusters\n",
    "\n",
    "    Returns:\n",
    "        purity:   a float, the purity of clustering result\n",
    "    \"\"\"\n",
    "    purity = 0.0\n",
    "    #YOUR CODE HERE\n",
    "    #YOUR CODE HERE\n",
    "    return purity\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eps_lst = np.linspace(0.1, 1, 19) #[0.1,0.15,...,1] for eps\n",
    "minpts_lst = np.linspace(2, 15, 14) #[2,3,...,15] for minpts\n",
    "purity_table = np.zeros((len(eps_lst), len(minpts_lst)))\n",
    "\n",
    "for i in range(len(eps_lst)):\n",
    "    for j in range(len(minpts_lst)):\n",
    "        cl = dbscan(X,eps_lst[i], minpts_lst[j]) #use your dbscan\n",
    "        ctable = contingency_table(cl, ground_truth) #use your contingency table function above\n",
    "        prt = purity(ctable, n) #use your purity function above\n",
    "        purity_table[i, j] = prt\n",
    "print(purity_table) #purity with different parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 6: DBSCAN with Gaussian kernel\n",
    "\n",
    "In this exercise, we are going to alter the DBSCAN algorithm to use another kernel.\n",
    "\n",
    "As stated in [Zaki, p.388], DBSCAN is a special case of DENCLUE.\n",
    "In particular, if we let $h=\\epsilon$ and $\\xi = minPts$ with a discrete kernel, then the two algorithms will yield the same result.\n",
    "\n",
    "We are now going to go a step in the other direction.\n",
    "Thas is, we will add a Gaussian kernel to the DBSCAN algorithm.\n",
    "The Gaussian kernel is defined as in Equation (1).\n",
    "\n",
    "$$\n",
    "K(\\mathbf{z})=\\frac{1}{(2 \\pi)^{d / 2}} \\exp \\left\\{-\\frac{\\mathbf{z}^{T} \\mathbf{z}}{2}\\right\\}\n",
    "\\qquad \\qquad (1)\n",
    "$$\n",
    "\n",
    "The implications in terms of the algorithm are the following:\n",
    "\n",
    "1. When selecting core points, they are now going to depend on a threshold $\\xi$ of the density estimates $\\hat f(x)$.\n",
    "2. The threshold $\\epsilon$ is now going to be compared against $K(\\frac{x-x_i}{h})$ when calculating neighborhoods and density connectedness.\n",
    "3. The value $h$ which was previously fixed to $\\epsilon$ is now going to be a parameter to the model.\n",
    "\n",
    "As before, the code below serves as inspiration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def K(z):\n",
    "    return  np.exp(-np.dot(z, z)/2) / np.sqrt(2*np.pi)\n",
    "\n",
    "def densityConnected(x, k, ...): # You might need more parameters here.\n",
    "    for neighbor in x.Nx:\n",
    "        # Connect (potentially recursively)\n",
    "        pass\n",
    "\n",
    "def gaussian_dbscan(X, e, xi, h):\n",
    "    \"\"\"\n",
    "        Args:\n",
    "            X:       Matrix of shape [n, d] with data points on the rows.\n",
    "            e:       Epsilon distance for neighborhood calculations.\n",
    "            m:       Minimum neighbors in epsilon neighborhood for a point to be a core point.\n",
    "            dist_fn: Distance function to be used.\n",
    "        \n",
    "        Returns:\n",
    "            clusters:   A vector of shape [n,] with integers, indicating cluster assignments.\n",
    "                        Let clusters[i] == -1 if point x_i is an outlier and a non-negative \n",
    "                        integer corresponding to the cluster index of point x_i otherwise.\n",
    "    \"\"\"\n",
    "    \n",
    "    # TODO, code here.\n",
    "    \n",
    "    core = []\n",
    "    for x in X:\n",
    "        # Compute neighborhoods and identify cores\n",
    "        pass\n",
    "\n",
    "    k = 0\n",
    "    for c in core:\n",
    "        # Assign unassigned cores to clusters\n",
    "        # using densityConnect\n",
    "        pass\n",
    "\n",
    "    clusters = []\n",
    "    for i in range(k):\n",
    "        # Build cluster indicator vector\n",
    "        pass\n",
    "    \n",
    "    # End TODO\n",
    "    \n",
    "    return clusters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the code below to run your algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test clustering on iris_PC dataset\n",
    "np.random.seed(0)\n",
    "\n",
    "slow = True\n",
    "if slow:\n",
    "    X, _ = load_t7()\n",
    "    n, _ = X.shape\n",
    "    X = X[np.random.permutation(n)]\n",
    "    X = X[:2000]\n",
    "    clustering = gaussian_dbscan(X, e=0.0001, xi=0.4, h=3.25)\n",
    "else: \n",
    "    X, _ = load_iris_PC()\n",
    "    clustering = gaussian_dbscan(X, e=0.2, xi=12., h=0.7)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "sc = ax.scatter(*(X[clustering >= 0].T), c=clustering[clustering>=0])\n",
    "sc = ax.scatter(*(X[clustering < 0].T), marker='x')\n",
    "plt.colorbar(sc)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 7: Comparing the DBSCAN versions above\n",
    "Please describe how the two versions of DBSCAN above differ.\n",
    "How do you think the the differences affects the clustering of the data?\n",
    "\n",
    "Did you see any practical differences in the \"experiments\" above?"
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}