05-clustering-hdbscan.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5-Clustering with HDBSCAN\n",
    "\n",
    "In this notebook we calculate the areas where the trips most commonly start and end. We call these the _trip endpoints_. Usually you also find these named as \"_stay points_\" in the literature but, due to the nature of the data that we are handling, where vehicles do not actuaally stop at these locations, we will keep the name \"_endpoint_\".\n",
    "\n",
    "**Requirements:**\n",
    "\n",
    "- Please run the `02-import-data.ipynb` notebook first, in case you need to buid the supporting SQLite database.\n",
    "- Required install: [ipywidgets](https://ipywidgets.readthedocs.io/en/stable/user_install.html). Enable using `jupyter nbextension enable --py widgetsnbextension --sys-prefix` for Jupyter Notebook and `jupyter labextension install @jupyter-widgets/jupyterlab-manager` for Jupyter Lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import folium\n",
    "import hdbscan\n",
    "import utm\n",
    "import ipywidgets as widgets\n",
    "\n",
    "from db.api import VedDb\n",
    "from h3 import h3\n",
    "from tqdm.notebook import tqdm\n",
    "from sklearn.cluster import DBSCAN\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "from folium.vector_layers import CircleMarker\n",
    "from colour import Color\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by instantiating the `VedDb` object to interface with the SQLite database created and prepared in previous notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "db = VedDb()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we collect the start and end locations of all the moves in the database as latitude and longitude pairs. Note that we have to join the `move` table to the `signal` table twice in order to get both start and end loctions. The result is downloaded into a DataFrame for in-memory manipulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sql = \"\"\"\n",
    "select     m.move_id\n",
    ",          si.latitude  as lat_ini\n",
    ",          si.longitude as lon_ini\n",
    ",          se.latitude  as lat_end\n",
    ",          se.longitude as lon_end\n",
    "from       move m\n",
    "inner join signal si on si.day_num = m.day_num and \n",
    "                        si.vehicle_id = m.vehicle_id and \n",
    "                        si.time_stamp = m.ts_ini\n",
    "inner join signal se on se.day_num = m.day_num and \n",
    "                        se.vehicle_id = m.vehicle_id and \n",
    "                        se.time_stamp = m.ts_end\n",
    "\"\"\"\n",
    "df_pt = db.query_df(sql)\n",
    "df_pt.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, for each move we collect the start location `(lat_ini, lon_ini)` and destination location `(lat_end, lon_end)`. Using this data, we can determine the implied geographical clusters using [HDBSCAN](https://hdbscan.readthedocs.io/en/latest/). Before running that algorithm, we must first collect all locations into a single array, keeping the original ordering. The first half of the array contains the start locations while the second part contains the destination locations. This is important in order to guarantee that we can assign the detetcted cluster identifiers back to the `move` table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc_ini = df_pt[['lat_ini', 'lon_ini']].to_numpy()\n",
    "loc_end = df_pt[['lat_end', 'lon_end']].to_numpy()\n",
    "locations = np.vstack((loc_ini, loc_end))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now cluster all the points. Here I am providing two functions to fit an HDBSCAN clusterer object. The `fit_utm_clusterer` function fits an HDBSCAN model using [UTM coordinates](https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_utm_clusterer(locations, min_cluster_size=20, min_samples=20):\n",
    "    xyzz = [utm.from_latlon(ll[0], ll[1]) for ll in locations]\n",
    "    pts = [[p[0], p[1]] for p in xyzz]\n",
    "\n",
    "    clusterer = hdbscan.HDBSCAN(min_cluster_size=min_cluster_size, \n",
    "                                min_samples=min_samples,\n",
    "                                metric='euclidean')\n",
    "    clusterer.fit(pts)\n",
    "    return clusterer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `fit_latlon_clusterer` uses the location coordinates directly and calculates distances using the haversine function. Both functions should produce equivalent results, but this one works better for larger areas, especially for longitude renges of over 6 degrees. For a small are such as is the case, the UTM version is adequate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_latlon_clusterer(locations, min_cluster_size=20, min_samples=20):\n",
    "    pts = np.radians(locations)\n",
    "\n",
    "    clusterer = hdbscan.HDBSCAN(min_cluster_size=min_cluster_size, \n",
    "                                min_samples=min_samples,\n",
    "                                metric='haversine')\n",
    "    clusterer.fit(pts)\n",
    "    return clusterer.labels_, clusterer.outlier_scores_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we run the code against the location array and collect the cluster identifiers in the same order as that of the inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusterer = fit_utm_clusterer(locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The list of unique cluster identifiers is calculated below by removing the first unique value (-1), the noise indicator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_clusters = np.unique(clusterer.labels_)[1:]\n",
    "print(\"The initial number of clusters is: {0}\".format(unique_clusters.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Detecting Outliers\n",
    "We are not done yet. The HDBSCAN algorithm can generate some odd-looking clusters with points that are seemingly misplaced, too far away from the main cluster. You will be able to judge that by yourself below while browsing through the map. Fortunately, HDBSCAN also provides a means for us to weed out such points using the concept of [outliers](https://hdbscan.readthedocs.io/en/latest/outlier_detection.html). For each point, the algorithm calculates an outlier score, where higher values mean higher likelihood of the point being an outlier.\n",
    "\n",
    "Outlier scores range from zero to one, and we can use such information to filter out the most glaring outliers. In this case, through visual inspection, I decided to filter out points with outllier scores above 0.8. Below, I try to provide a more quantitative support for this decision.\n",
    "\n",
    "Let's start by plotting a histogram of the outlier score distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hh = plt.hist(clusterer.outlier_scores_, bins=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, there are very few outlying points beyond _0.8_. The proportion of points with outlier scores below _0.8_ is quite large indeed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores = clusterer.outlier_scores_\n",
    "print(\"Below 0.8 you can find {:.2f}% of the cluster points.\".format(100.0 * scores[scores < 0.8].shape[0] / scores.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have another look at this distribution. In the following plot, the _x_ axis reepresents the outlier score and the _y_ axis represents the proportion of cluster points with up to the given outlier score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [8, 8]\n",
    "score_shape = [scores[scores < x].shape[0] / scores.shape[0] for x in np.arange(0.0, 1.0, 0.01)]\n",
    "plt.plot(np.arange(0.0, 1.0, 0.01), score_shape)\n",
    "plt.xlabel(\"Outlier Score\")\n",
    "plt.ylabel(\"Point Proportion\")\n",
    "plt.title(\"Outlier Score Cumulative Distribution\")\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.Series(scores).quantile(0.95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see some visuals!\n",
    "\n",
    "## Interactive Cluster Exploration\n",
    "\n",
    "Here you can explore all the generated clusters through the interactive widget below. The `show_cluster_map` displays a cluster map with all the points colored according to their outlier score. The more a point color shifts to red, the higher its outlier score. For your convenience, each point has a tooltip with its score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_cluster_map(cluster_id):\n",
    "    blue = Color(\"blue\")\n",
    "    red = Color(\"red\")\n",
    "    color_range = list(blue.range_to(red, 10))\n",
    "    \n",
    "    map = folium.Map(prefer_canvas=True, tiles='CartoDB positron')\n",
    "    \n",
    "    clusters = clusterer.labels_\n",
    "    outlier_scores = clusterer.outlier_scores_\n",
    "    \n",
    "    points = locations[clusters == cluster_id]\n",
    "    scores = outlier_scores[clusters == cluster_id]\n",
    "    \n",
    "    for i in range(points.shape[0]):\n",
    "        point = points[i]\n",
    "        color = color_range[int(scores[i] * 10)]\n",
    "        CircleMarker(point, radius=1, color=color.hex, tooltip=\"{:.2f}\".format(scores[i])).add_to(map)\n",
    "        \n",
    "    min_lat, max_lat = points[:, 0].min(), points[:, 0].max()\n",
    "    min_lon, max_lon = points[:, 1].min(), points[:, 1].max()\n",
    "    map.fit_bounds([[min_lat, min_lon], [max_lat, max_lon]])\n",
    "    return map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ii = interact(show_cluster_map, cluster_id=widgets.IntSlider(min=0, max=clusterer.labels_.max(), step=1, value=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outlier Filtering\n",
    "\n",
    "We now must make good on our promise and filter the offending outliers. We do this by marking as noise all the points with an outlier score above _0.8_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_clusters = np.copy(clusterer.labels_)\n",
    "filtered_clusters[scores >= 0.8] = -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Serialize to the Database\n",
    "\n",
    "For the sake of future convenience, we will now create and fill in a table whith all the cluster points. This table will make our life much easier in the future, when handling clusters, their locations and shapes.\n",
    "\n",
    "We start by making sure that the table exists and is empty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not db.table_exists(\"cluster_point\"):\n",
    "    sql = \"\"\"\n",
    "    CREATE TABLE cluster_point (\n",
    "        pt_id           INTEGER PRIMARY KEY ASC,\n",
    "        cluster_id      INT NOT NULL,\n",
    "        latitude        FLOAT NOT NULL,\n",
    "        longitude       FLOAT NOT NULL,\n",
    "        h3              TEXT\n",
    "    )\n",
    "    \"\"\"\n",
    "    db.execute_sql(sql)\n",
    "    db.execute_sql(\"CREATE INDEX idx_cluster_point_cluster ON cluster_point (cluster_id)\")\n",
    "else:\n",
    "    db.execute_sql(\"DELETE FROM cluster_point\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can insert the locations for each cluster, along with the [H3](https://eng.uber.com/h3/) hexagon codes at [resolution level 12](https://uber.github.io/h3/#/documentation/core-library/resolution-table). In the next notebook, we will display the clusters using the outline of the collated shapes of all cluster hexagons.\n",
    "\n",
    "To prepare the insert statement, we now collect all the input data into a list:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h3_level = 12\n",
    "cluster_points = []\n",
    "for i in tqdm(range(filtered_clusters.shape[0])):\n",
    "    if filtered_clusters[i] >= 0:\n",
    "        pt = (int(filtered_clusters[i]), locations[i, 0], locations[i, 1], h3.geo_to_h3(locations[i, 0], locations[i, 1], h3_level))\n",
    "        cluster_points.append(pt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now insert the cluster points into the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "db.insert_cluster_points(cluster_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DataFrame can also be updated with the cluster identifiers like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = filtered_clusters.shape[0] // 2\n",
    "df_pt['cluster_ini'] = filtered_clusters[:n]\n",
    "df_pt['cluster_end'] = filtered_clusters[n:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how the DataFrame looks like after cluster identifier assignment. Note that the constant `-1` means that the point was not assigned to any cluster, and was considered as noise instead. We will exclude these from the future trajectory analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pt.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying the `move` Table\n",
    "\n",
    "Now that we have the clusters identified, we can assign them back to the `move` table. Unfortunately, this table has no columns where to store the clusters identifiers, so we must first handle that. For convenience, we also create an extra index on the new columns for more convenient search on them.\n",
    "\n",
    "**Notes**: \n",
    "- The function `table_has_column` tests whether a given table has a named column. Here we only test against one column as it is enough.\n",
    "- Depending on the time you use this code, the database may already have been created with these columns and indexes. If that is the case, the code below does nothing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not db.table_has_column('move', 'cluster_ini'):\n",
    "    db.execute_sql(\"alter table move add cluster_ini INT not null default -1\")\n",
    "    db.execute_sql(\"alter table move add cluster_end INT not null default -1\")\n",
    "    db.execute_sql(\"create index idx_move_clusters on move (cluster_ini, cluster_end)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to update the `move` table using the recently calculated cluster identifiers. To do so, we use the data stored in the DataFrame to feed an update query. Note that we retrieve the data as a list of tuples with the eaxct order for consumption in the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clids = list(df_pt[['cluster_ini', 'cluster_end', 'move_id']].itertuples(index=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "db.update_move_clusters(clids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now associated each move with a pair of endpoint clusters, and this allows us to perform more powerful analyses to the data, such as determine how many trips occur between two general endpoints, how many different trajectories there are and what vehicles have used them and at what time, consuming how much fuel or energy.\n",
    "\n",
    "Now, we turn to the issue of generating a geo-fence for the clusters so they can be easily displayed on a map."
   ]
  },
  {
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}