Jupyter/M06-Advanced/M06B-IBMDP-Exploration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Cloud-First](../image/CloudFirst.png) \n",
    "\n",
    "# SIT742: Modern Data Science \n",
    "**(Module: Advanced Topics in Big Data)**\n",
    "\n",
    "---\n",
    "- Materials in this module include resources collected from various open-source online repositories.\n",
    "- You are free to use, change and distribute this package.\n",
    "- If you found any issue/bug for this document, please submit an issue at [tulip-lab/sit742](https://github.com/tulip-lab/sit742/issues)\n",
    "\n",
    "\n",
    "Prepared by **SIT742 Teaching Team**\n",
    "\n",
    "---\n",
    "\n",
    "\n",
    "## Session 6B: IBM DP Library for Exploratory Data Analysis\n",
    "\n",
    "This set of notebooks is from [IBM Diffprivlib](https://github.com/IBM/differential-privacy-library/tree/main/notebooks), and you may want to directly access those other sample notebooks for its usage.\n",
    "\n",
    "**The purpose of this session is to illustrate**\n",
    "\n",
    "* How to use IBM DP Library for data exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import diffprivlib as dp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data and differential privacy set-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by initialising a `BudgetAccountant` to track the privacy spend of our exploration. By using the `set_default()` method, this accountant will be used by default in all diffprivlib calls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = dp.BudgetAccountant(1, 0)\n",
    "acc.set_default()\n",
    "\n",
    "# Our default epsilon value for this exploration\n",
    "eps = 0.04"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the Covertype dataset for this example. We can load is using `sklearn`. For the purpose of this exploration, we will use a numerical encoding of the \"Wilderness Area\" and \"Soil Type\" attributes. This will simplify visualisations later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import fetch_covtype\n",
    "\n",
    "dataset = fetch_covtype()\n",
    "data, labels = dataset.data, dataset.target\n",
    "\n",
    "# Switch from one-hot encoding to numerical encoding for Wilderness Area and Soil Type\n",
    "data = np.concatenate((\n",
    "    data[:, :-44], \n",
    "    (data[:, -44:-40] * np.arange(0, 4)).sum(axis=1)[:, np.newaxis],\n",
    "    (data[:, -40:] * np.arange(0, 40)).sum(axis=1)[:, np.newaxis]\n",
    "), axis=1)\n",
    "\n",
    "n_examples = data.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To assist our exploration, we will also specify the column names, and the ranges of each column. In this case, the ranges are assumed to be public knowledge, and therefore are not required to be calculated on the data. If the range is not known, it is common to search using a differentially private histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col_names = [\"Elevation\", \"Aspect\", \"Slope\", \"H-Distance To Hydrology\", \n",
    "             \"V-Distance To Hydrology\", \"H-Distance To Roadways\", \n",
    "             \"Hillshade 9am\", \"Hillshade Noon\", \"Hillshade 3pm\", \n",
    "             \"H-Distance To Fire Points\", \"Wilderness Area\", \"Soil Type\"]\n",
    "\n",
    "ranges = [(1000.0, 4000.0), (0, 360), (0, 66), (0, 1000.0), (0, 1000.0),\n",
    "  (10, 10000.0), (0, 255), (0, 255), (0, 255), (10, 10000.0), (0, 3), (0, 39)]\n",
    "\n",
    "assert len(col_names) == len(ranges)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `histogram` function is a useful tool to visualise the data in a differentially private way. The syntax is the same as the corresponding function in NumPy, with the addition of an `epsilon` parameter. Also, to preserve differential privacy, we must specify the `range` of the histogram.\n",
    "\n",
    "In this first example, we plot the number of examples of each of the Covertype labels, 1‒7. The counts are not exact due to differential privacy, but give a good approximation of the true values. We can see that almost 50% of examples in the dataset are associated with the label \"2\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(5,3.5))\n",
    "\n",
    "ax.bar(np.linspace(1, 7, 7), dp.tools.histogram(labels, epsilon=eps, bins=7, range=(1, 8))[0] / n_examples)\n",
    "ax.set_title(\"Label counts\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the distribution of features in the dataset using `histogram`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(3, 4, sharey=True, figsize=(16, 12))\n",
    "fig.subplots_adjust(hspace=0.35)\n",
    "\n",
    "for i in range(len(col_names)):\n",
    "    ax = axs.flat[i]\n",
    "    ax.set_title(col_names[i])\n",
    "    _range = ranges[i]\n",
    "    n_bins = 10\n",
    "    \n",
    "    # Wilderness Area and Soil Type need special config\n",
    "    if i >= 10:\n",
    "        _range = (0, _range[1] + 1)\n",
    "        n_bins = _range[1]\n",
    "\n",
    "    bar_width = (_range[1] - _range[0]) / n_bins * 0.9\n",
    "    \n",
    "    hist, bins = dp.tools.histogram(data[:, i], bins=n_bins, range=_range, epsilon=eps)\n",
    "    ax.bar(bins[:-1] + np.diff(bins), hist / n_examples, width=bar_width)\n",
    "    ax.set_ylim(0, 1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After our initial exploration, now is a good time to examine our budget spend, confirming that we have executed 13 queries thus far."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "acc.total()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(acc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use 2-dimensional histograms to extract more information in a differentially private way. In this example, we can plot the distribution of a feature with respect to each label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(15, 4), sharey=True)\n",
    "n_bins = 10\n",
    "\n",
    "# Columns of the dataset we want to plot\n",
    "cols = (1, 3, 8)\n",
    "\n",
    "for _i, col in enumerate(cols):\n",
    "    ax = axs.flat[_i]\n",
    "    ax.set_title(col_names[col])\n",
    "    width = (ranges[col][1] - ranges[col][0]) / n_bins * 0.9\n",
    "\n",
    "    hist, bins, _ = dp.tools.histogram2d(data[:, col], labels, epsilon=eps, bins=(n_bins, 7), \n",
    "                                             range=(ranges[col], (1, 8)))\n",
    "    cum_sum = np.zeros_like(hist[:, 0])\n",
    "\n",
    "    for i in range(hist.shape[1]):\n",
    "        ax.bar(bins[:-1] + np.diff(bins), hist[:, i] / n_examples, bottom=cum_sum, label=str(i+1), width=width)\n",
    "        cum_sum += hist[:, i] / n_examples\n",
    "        \n",
    "ax.legend(title=\"Cover type\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With some extra effort, we can also produce 2-dimensional histograms to show correlation between features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(5, 4))\n",
    "col1, col2 = (3, 5)\n",
    "\n",
    "n_bins = 10\n",
    "\n",
    "ax.set_title(col_names[col1] + \" v \\n\" + col_names[col2])\n",
    "width = (ranges[col1][1] - ranges[col1][0]) / n_bins * 0.9\n",
    "\n",
    "hist, bins, _ = dp.tools.histogram2d(\n",
    "    data[:, col1], \n",
    "    (data[:, col2] > 1000).astype(int) + (data[:, col2] > 2000) + (data[:, col2] > 3000) + (data[:, col2] > 4000), \n",
    "    epsilon=eps, bins=(n_bins, 5), range=(ranges[col1], (0, 5))\n",
    ")\n",
    "legends = [\"<= 1000\", \"> 1000 && <= 2000\", \"> 2000 && <= 3000\", \"> 3000 && <= 4000\", \"> 4000\"]\n",
    "\n",
    "cum_sum = np.zeros_like(hist[:, 0])\n",
    "\n",
    "for i in range(hist.shape[1]):\n",
    "    ax.bar(bins[:-1] + np.diff(bins), hist[:, i] / n_examples, bottom=cum_sum, label=legends[i], width=width)\n",
    "    cum_sum += hist[:, i] / n_examples\n",
    "\n",
    "ax.set_xlabel(col_names[col1])\n",
    "ax.set_ylabel(\"Frequency\")\n",
    "ax.legend(title=col_names[col2])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Colour maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the 2-dimensional histogram to plot colour maps of the data. This allows us to examine correlations between features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "col1, col2 = (1, 3)\n",
    "\n",
    "hist, bins_y, bins_x = dp.tools.histogram2d(data[:, col1], data[:, col2], epsilon=eps, \n",
    "                                            range=(ranges[col1], ranges[col2]))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(5, 4))\n",
    "c = ax.pcolormesh(bins_y, bins_x, hist / n_examples, cmap='Blues')\n",
    "ax.set_ylabel(col_names[col1])\n",
    "ax.set_xlabel(col_names[col2])\n",
    "fig.colorbar(c, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other queries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Diffprivlib includes a number of simple statistical functions to examine the data in more detail. These include `mean`, `var`, `sum` and `count_nonzero` functions.  They can be useful in querying more specific information than can be portrayed by histograms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dp.tools.mean(data[:, 0], bounds=ranges[0], epsilon=eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp.tools.var(data[:, 0], bounds=ranges[0], epsilon=eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`count_nonzero` is useful for counting more specific information than is required from a histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp.tools.count_nonzero(data[:, 0] >= 3000, epsilon=eps) / n_examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get label-specific information by indexing the `data` array appropriately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp.tools.mean(data[labels==1, 0], bounds=ranges[0], epsilon=eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total privacy loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having completed our initial exploration, we can compute the total privacy loss from the accountant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc.total()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can also compute the remaining budget to be spent on subsequent queries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acc.remaining()"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}