smd_unsupervised.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:26.174916Z",
     "start_time": "2018-11-27T15:05:26.151988Z"
    }
   },
   "outputs": [],
   "source": [
    "%%javascript\n",
    "$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unsupervised Machine Learning  and Clustering\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2 id=\"tocheading\">Table of Contents</h2>\n",
    "<div id=\"toc\"></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The story so far:\n",
    "\n",
    "- Linear Discriminant Analysis (LDA) and Fisher's linear discriminant\n",
    "- Principal Component Analysis (PCA)\n",
    "- Feature Selection\n",
    "\n",
    "\n",
    "Notations:\n",
    "\n",
    "* Uppercase letters such as $X$ or $Y$ denote generic aspects of a variable (i.e. the actual random variable)\n",
    "* Observed values are written in lowercase. The ith observed value of $X$ is written as $x_i$\n",
    "* Matrices are written in bold uppercase letters as in $\\mathbf{X}$\n",
    "* Observations map as *rows* in the matrix while the observed variables are the *columns*.\n",
    "\n",
    "So if I measure two observables $p = 2$ the size and weight of $N = 100$ people, I get a $N \\times p$ matrix $\\mathbf{X}$.\n",
    "One observation in that matrix is denoted as $x_i = [ size, weight ]$ while all observations of the variable size are denoted by $\\mathbf{x}_j$ \n",
    "\n",
    "\n",
    "So far we have been occupied with \n",
    "predicting the values of one or more outputs or response variables $Y = (Y_1, \\ldots, Y_m)$ for a given set of input or predictor variables $X = (X_1, \\ldots , X_p)$. \n",
    "\n",
    "One possible definition of supervised machine learning\n",
    "\n",
    "> Given a $N \\times p$ matrix $\\mathbf{X}$ and some associated output vector $\\mathbf{Y} \\in \\mathbb{R}^N$,\n",
    " find a function $f(X) = \\hat{Y}$ that takes a vector $X \\in \\mathbb{R}^p$ and returns a prediction for $\\hat{Y}$\n",
    " where some \"loss function\" $L(Y, f(X))$ is minimized for all $X$.\n",
    " \n",
    "We've tried to split points into different classes by finding a decision function which seperates the points in an \"optimal\" way. This process is often called __classification__\n",
    "\n",
    "In __regression__ the dependent variable is not discrete but a continous value. The problem remains similar however. From *known* input data we try to find a function which accurately predicts $Y \\in \\mathbb{R}$.\n",
    "\n",
    "In __unsupervised__ machine learning problems the output vector $\\mathbf{y}$ is *unknown*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:27.643365Z",
     "start_time": "2018-11-27T15:05:26.177785Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from ml import plots\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib.colors import ListedColormap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:27.643365Z",
     "start_time": "2018-11-27T15:05:26.177785Z"
    }
   },
   "outputs": [],
   "source": [
    "pd.options.display.max_rows = 10\n",
    "plots.set_plot_style()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification Example (Iris)\n",
    "\n",
    "The sterotypical machine learning data set. Given some measurements from flower petals is it possible to deduce the species for a new sample of measurements?\n",
    "\n",
    "Originally published in 1936 by Ronald Fisher in his paper about LDA. \n",
    "\n",
    "<table style=\"width:100%\">\n",
    "  <tr>\n",
    "    <td><img src=\"./ml/images/setosa.jpg\" alt=\"Drawing\" style=\"width: 200px;\"/></td>\n",
    "    <td><img src=\"./ml/images/versicolor.jpg\" alt=\"Drawing\" style=\"width: 200px;\"/></td> \n",
    "    <td><img src=\"./ml/images/virginica.jpg\" alt=\"Drawing\" style=\"width: 200px;\"/></td>\n",
    "  </tr>\n",
    "</table>\n",
    "\n",
    "<p style=\"color:gray\"> Iris flower data set. (2017, October 24). In Wikipedia, The Free Encyclopedia.  <p>\n",
    "\n",
    "Length and width of the petals and sepals making the dataset four dimensional.  Measurements where taken in Canada by Botanist Edgar Anderson:\n",
    "\n",
    "> all from the same pasture, and picked on the same day and measured at the same time by the same person with the same apparatus\n",
    "\n",
    "The label vector $\\mathbf{y}$ is known because the botanist categorized the flowers into species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:28.673241Z",
     "start_time": "2018-11-27T15:05:27.646648Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_iris\n",
    "from itertools import combinations\n",
    "\n",
    "\n",
    "\n",
    "iris = load_iris()\n",
    "X = iris.data  # we only take the first two features.\n",
    "y = iris.target\n",
    "\n",
    "\n",
    "\n",
    "discrete_cmap = ListedColormap([f'C{i}' for i in range(len(iris.target_names))])\n",
    "\n",
    "f, axs = plt.subplots(3, 2, constrained_layout=True)\n",
    "axs = axs.flatten()\n",
    "\n",
    "for i, (col_x, col_y) in enumerate(combinations(range(X.shape[1]), 2)):\n",
    "\n",
    "    plot = axs[i].scatter(\n",
    "        X[:, col_x], X[:, col_y],\n",
    "        c=y,\n",
    "        cmap=discrete_cmap,\n",
    "        label='species',\n",
    "        vmin=-0.5, vmax=len(iris.target_names) - 0.5,\n",
    "    )\n",
    "    axs[i].set_xlabel(iris.feature_names[col_x].replace('(cm)', '/ cm'))\n",
    "    axs[i].set_ylabel(iris.feature_names[col_y].replace('(cm)', '/ cm'))\n",
    "    \n",
    "bar = f.colorbar(plot)\n",
    "bar.set_ticks(np.arange(len(iris.target_names)))\n",
    "bar.set_ticklabels(iris.target_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The question this data set can help us answer is:\n",
    "\n",
    "> What species does this flower likely belong to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised Learning\n",
    "\n",
    "In unsupervised learning we have no given $Y$.\n",
    "\n",
    "These methods try to find the underlying (joint) probability density $Pr(X)$ so that we might learn some properties about it.\n",
    "\n",
    "One common question is whether $X$ is created by a mixture of two or more underlying random variables.\n",
    "\n",
    "In simpler terms:\n",
    "\n",
    " - Given some points, can we infer something about the underlying distributions or labels\n",
    "\n",
    "\n",
    "### Clustering\n",
    "\n",
    "Clustering algorithms try to find modes of $Pr(X)$ based on densities, neighbourhood relations or any other measure of 'similarity'  between points.\n",
    "\n",
    "More generally speaking to quote wikipedia again:\n",
    "\n",
    ">Cluster analysis or clustering is the task of grouping a set of objects in such a way that objects in the same group (called a cluster) are more similar (in some sense or another) to each other than to those in other groups (clusters).\n",
    "\n",
    "\n",
    "In the example below an unknown number of gaussian distributions are sampled to create a scatter plot.\n",
    "\n",
    "Can one infer $Pr(X)$ by looking at $X$, i.e. the blue dots?\n",
    "\n",
    "In this case we know that this distribution of blue dots, $X$, was created by sampling $k$ two dimensional gaussians with known standard deviation.\n",
    "This is just what \n",
    "\n",
    "    X, y = make_blobs(n_samples=300, centers=k, center_box=(-2, 2), cluster_std=0.5)\n",
    " \n",
    "does.\n",
    "We even know in what region of space we have to look for to find the centroids centroids of these blobs. \n",
    "\n",
    "Questions: \n",
    "\n",
    "   - How many blobs do we have?\n",
    "   - Which points belong to which?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:28.984232Z",
     "start_time": "2018-11-27T15:05:28.675757Z"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(1234)\n",
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "# choose a random number between 1 and 3\n",
    "k = np.random.randint(1, 4)\n",
    "\n",
    "# create k blobs \n",
    "X, y = make_blobs(n_samples=300, centers=k, center_box=(-2, 2), cluster_std=0.5)\n",
    "\n",
    "# plot points without any color coding.\n",
    "plt.scatter(X[:, 0], X[:, 1])\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "If you guessed k=2, you're wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:29.586188Z",
     "start_time": "2018-11-27T15:05:28.987673Z"
    }
   },
   "outputs": [],
   "source": [
    "# plot the same points this time with color coding \n",
    "plt.scatter(X[:, 0], X[:, 1], c=y, cmap=discrete_cmap)\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### K-Means Algorithm\n",
    "\n",
    "The k-Means algorithm tries to find a partition of the data into $k$ clusters $S = \\{S_1, \\ldots, S_k\\}$  which minimize the variance within those clusters. The number of clusters $k$ has to specified by the user.\n",
    "\n",
    "Formally speaking the algorithm solves\n",
    "$$\n",
    "{\\underset {S}{\\operatorname {arg\\,min} }}\\sum _{S_i \\in S}\\sum _{x \\in S_{i}}\\left\\|x -{\\overline{x}}_{S_i}\\right\\|^{2}.\n",
    "$$\n",
    "\n",
    "\n",
    "This problem is similar in nature to that of feature selection. A minimization over all possible subsets of the points $S$. Solving the general optimization problem is NP-Hard and therefore intractable. \n",
    "\n",
    "In simpler terms, it would take too long to find the optimal solution in the general case. \n",
    "Again there is a popular greedy heuristic which is usualy used to solve the problem.\n",
    "\n",
    "#### Loyds Algorithm\n",
    "\n",
    "Loyds Algorithm (sometimes also simply called *the* k-means algorithm) finds a local optimum using a greedy heuristic.\n",
    "\n",
    "It does so iterativly according to the following steps \n",
    "\n",
    "1. Pick some initial cluster means (or centroids) $\\{m_1, \\ldots, m_k \\}$ either randomly or according to some heuristic.\n",
    "\n",
    "2. Create a partition $S$  by assigning each point $x \\in X$ to the cluster $S_i$ where the distance to $m_i$ is the smallest.\n",
    "\n",
    "3. Update the cluster means by calculating the means within the assigned clusters. \n",
    "\n",
    "4. Repeat steps 2 and 3 until convergence.\n",
    "\n",
    "One popular convergence is the total change in cluster centroids per iteration. If the cluster centroids do not change from one iteration to the next, return the cluster centroids.\n",
    "\n",
    "\n",
    "The example below shows a k-mean clustering on $k=3$ random blobs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:30.581347Z",
     "start_time": "2018-11-27T15:05:29.591438Z"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(1234)\n",
    "from sklearn.cluster import KMeans\n",
    "\n",
    "# create three random blobs\n",
    "X, y = make_blobs(n_samples=300, centers=3, center_box=(-2, 2), cluster_std=0.5)\n",
    "\n",
    "# use KMeans to predict clusters for each sample in X\n",
    "prediction = KMeans(n_clusters=3).fit_predict(X)\n",
    "\n",
    "# shift labels to get the right colors\n",
    "prediction = np.choose(prediction, [2, 0, 1])\n",
    "\n",
    "# plot shades with predicted clusters\n",
    "plt.scatter(X[:, 0], X[:, 1], c=prediction, cmap=discrete_cmap, alpha=0.3, s=250, lw=0, label='prediction')\n",
    "\n",
    "# plot points with true cluster associations\n",
    "plt.scatter(X[:, 0], X[:, 1], c=y, cmap=discrete_cmap, label='truth')\n",
    "plt.legend(loc='upper right')\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below you can see an animation of k-Means as it converges. The black hexagons indicate the cluster centroids in the corresponding iteration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:58.415194Z",
     "start_time": "2018-11-27T15:05:30.584504Z"
    }
   },
   "outputs": [],
   "source": [
    "from IPython import display\n",
    "import time\n",
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "np.random.seed(1)\n",
    "\n",
    "# create 4 clsuters\n",
    "k = 4\n",
    "n_iters = 25\n",
    "discrete_cmap = ListedColormap([f'C{i}' for i in range(k)])\n",
    "\n",
    "# choose inital cluster centers\n",
    "X, y = make_blobs(n_samples=500, centers=k, center_box=(-2, 2), cluster_std=.5, random_state=1234)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.set_axis_off()\n",
    "ax.set_xlim(-3, 3)\n",
    "ax.set_ylim(-3, 3)\n",
    "\n",
    "\n",
    "init_centers = np.random.uniform(-6, 6, size=[4, 2])\n",
    "center_history = np.zeros((n_iters, k, 2))\n",
    "\n",
    "\n",
    "center_lines = [ax.plot([], [])[0] for i in range(k)]\n",
    "points = ax.scatter(X[:,0], X[:,1], c='k', cmap=discrete_cmap, alpha=0.5)\n",
    "center_plot = ax.scatter(\n",
    "    init_centers[:, 0],\n",
    "    init_centers[:, 1],\n",
    "    c=np.arange(k),\n",
    "    cmap=discrete_cmap,\n",
    "    marker='h',\n",
    "    edgecolor='k',\n",
    "    s=400,\n",
    "    label='cluster center',\n",
    ")\n",
    "\n",
    "ax.legend(loc='upper right')   \n",
    "    \n",
    "    \n",
    "def update(i):\n",
    "    kmeans = KMeans(n_clusters=k, init=init_centers, max_iter=i + 1, n_init=1)\n",
    "    prediction = kmeans.fit_predict(X)\n",
    "    \n",
    "    center_history[i] = kmeans.cluster_centers_\n",
    "    \n",
    "    for j, line in enumerate(center_lines):\n",
    "        line.set_data(center_history[:i + 1, j, 0], center_history[:i + 1, j, 1])\n",
    "    \n",
    "    center_plot.set_offsets(kmeans.cluster_centers_)\n",
    "    points.set_array(prediction)\n",
    "    points.set_cmap(discrete_cmap)\n",
    "    ax.set_title('Iteration {}'.format(i + 1))\n",
    "    \n",
    "\n",
    "\n",
    "ani = FuncAnimation(fig, update, blit=False, frames=n_iters)\n",
    "plt.close()\n",
    "display.HTML(ani.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the number of clusters in the center is not known, the algorithm does not produce meaningful results.\n",
    "This is a serious limitation as the number of clusters in the data is rarely known.  \n",
    "\n",
    "\n",
    "In the example below we sample two uniform 2D distributions and use k-means to cluster them into $k=4$ regions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:58.745695Z",
     "start_time": "2018-11-27T15:05:58.417540Z"
    }
   },
   "outputs": [],
   "source": [
    "u1 = np.random.uniform(-1, 0, size=(500, 2))\n",
    "u2 = np.random.uniform(0, 1, size=(500, 2))\n",
    "X = np.append(u1, u2, axis=0)\n",
    "plt.scatter(X[:, 0], X[:, 1], color='gray')\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:59.247113Z",
     "start_time": "2018-11-27T15:05:58.748135Z"
    }
   },
   "outputs": [],
   "source": [
    "# use KMeans to predict clusters for each sample in X using the deliberatley wrong value of k=4\n",
    "prediction = KMeans(n_clusters=4).fit_predict(X)\n",
    "\n",
    "plt.scatter(X[:, 0], X[:, 1], c=prediction, cmap=discrete_cmap)\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:05:59.660863Z",
     "start_time": "2018-11-27T15:05:59.249232Z"
    }
   },
   "outputs": [],
   "source": [
    "# use KMeans to predict clusters for each sample in X this time using the correct value for k=2\n",
    "prediction = KMeans(n_clusters=2).fit_predict(X)\n",
    "\n",
    "plt.scatter(X[:, 0], X[:, 1], c=prediction, cmap=discrete_cmap)\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The k-Means algorithm works well on convex clusters with similar standard deviations. But it fails on elongated or concave shapes.\n",
    "\n",
    "Like all cluster algorithms k-Means has a number of advantages and disadvantages\n",
    "\n",
    "Pros:\n",
    "\n",
    " - Relatively fast\n",
    " - Predictable results when inital cluster centroids are fixed.\n",
    " - works well on convex and 'blob' like shapes\n",
    "\n",
    "Cons:\n",
    "\n",
    " - Number of clusters has to be known beforehand\n",
    " - Gets worse and slower in high dimensions\n",
    " \n",
    "The follwing snippet transforms two random blobs to be elongated along the x axis.\n",
    "In this case k-Means fails as the neighborhood relation does not yield meaningful results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:00.102111Z",
     "start_time": "2018-11-27T15:05:59.663065Z"
    }
   },
   "outputs": [],
   "source": [
    "X, y = make_blobs(n_samples=1300, centers=2, random_state=2)\n",
    "transformation = [[20, 0], [0, 1]]\n",
    "X = np.dot(X, transformation)\n",
    "prediction = KMeans(n_clusters=2,).fit_predict(X)\n",
    "\n",
    "plt.scatter(X[:, 0], X[:, 1], c=prediction, cmap=discrete_cmap)\n",
    "plt.gca().set_aspect('equal', 'datalim')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Scaling\n",
    "\n",
    "k-Means and many other algorithms, especially the ones using kernel functions, are sensitive to data scaling and centering.\n",
    "\n",
    "The scikit-learn toolkit includes a  preprocessing module with  methods to scale and transform data.\n",
    "\n",
    "##### Standard Scaling\n",
    "\n",
    "The simplest scaler is the `StandardScaler` scaling data to unit variance and zero mean. It makes it 'gaussian like'.\n",
    "\n",
    "##### MinMax Scaling\n",
    "\n",
    "Forces the values of each attribute in the data to be within the given feature range. When transforming  to values between 0 and 1 the applied transformation is\n",
    "\n",
    "$$\n",
    "\\mathbf{x}_j^{\\prime} = \\frac{(\\mathbf{x}_j - \\min(\\mathbf{x}_j))}{(\\max(\\mathbf{x}_j) - \\min(\\mathbf{x}_j))}\n",
    "$$\n",
    "\n",
    "\n",
    "##### Quantile Scaling\n",
    "\n",
    "Transform the given data to follow a normal, uniform or any other distribution with a CDF that can be inverted.\n",
    "Use inverse transform sampling to transform the data. I.e. apply the inverse CDF to each column in the data.\n",
    "\n",
    "https://en.wikipedia.org/wiki/Inverse_transform_sampling\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:01.555566Z",
     "start_time": "2018-11-27T15:06:00.104848Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn import preprocessing\n",
    "\n",
    "f, axs = plt.subplots(4, 1, figsize=(10, 20), sharex=True)\n",
    "\n",
    "scalers=[\n",
    "        preprocessing.StandardScaler(),\n",
    "        preprocessing.MinMaxScaler(feature_range=(-1, 1)),\n",
    "        preprocessing.MaxAbsScaler(),\n",
    "        preprocessing.QuantileTransformer()\n",
    "]\n",
    "\n",
    "for scaler, ax in zip(scalers, axs):\n",
    "    X_prime = scaler.fit_transform(X)\n",
    "    prediction = KMeans(n_clusters=2, random_state=0).fit_predict(X_prime)\n",
    "    ax.set_title(scaler.__class__.__name__)\n",
    "    ax.scatter(X_prime[:, 0], X_prime[:, 1], c=prediction, cmap=discrete_cmap)\n",
    "    \n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image Compression Using k-Means\n",
    "\n",
    "For a typical color, or RGB,  image the value for each color component must be stored. It is basically a list of color vectors \n",
    "\n",
    "$$\n",
    "c_i = \\begin{pmatrix} R, G, B \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "For a $420 \\times 640$ pixel image that results in  $420 \\cdot 640 \\cdot 3 \\cdot 2  \\; \\text{byte} = 1.6128 \\;  \\text{megabyte}$ of memory.\n",
    "This is how a typical Bitmap Image works.\n",
    "\n",
    "The GIF format allows the compression of color information.  \n",
    "\n",
    "In the example below we use k-Means to compress the color information of an image. We cluster the space of color components (either RGB, HSV or similar) into $k$ cluster centroids. \n",
    "Then only the association of each pixel to the cluster centroids and the centroids themselves have to be stored. \n",
    "\n",
    "When clustering the colors of an image into $k=50$ centroids for example  we only need to store $(420 \\cdot 640 + k) \\cdot 2  \\; \\text{byte} = 537.7 \\; \\text{kilobyte}$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:03.176783Z",
     "start_time": "2018-11-27T15:06:01.557800Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_sample_image\n",
    "from skimage import color\n",
    "\n",
    "#load an example image as a 3D array\n",
    "image_rgb = load_sample_image(\"flower.jpg\")\n",
    "image_rgb = np.array(image_rgb, dtype=np.float64) / 255\n",
    "\n",
    "#store width length and number of colors\n",
    "width, length, channels = image_rgb.shape  \n",
    "\n",
    "# show image\n",
    "plt.imshow(image_rgb)\n",
    "\n",
    "#convert image to hsv for nicer plots\n",
    "image_hsv = color.rgb2hsv(image_rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:28.343954Z",
     "start_time": "2018-11-27T15:06:03.179701Z"
    }
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits import mplot3d\n",
    "\n",
    "# plot H, V and S values of each pixel into a 3D coordinate system.\n",
    "\n",
    "fig = plt.figure(figsize=(12, 12))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "\n",
    "ax.scatter3D(\n",
    "    image_hsv[:, :, 0],\n",
    "    image_hsv[:, :, 1],\n",
    "    image_hsv[:, :, 2],\n",
    "    c=image_rgb.reshape(-1, 3),\n",
    "    alpha=0.5,\n",
    "    lw=0,\n",
    ")\n",
    "\n",
    "ax.set_xlabel('hue')\n",
    "ax.set_ylabel('saturation')\n",
    "ax.set_zlabel('lightness value')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:53.066282Z",
     "start_time": "2018-11-27T15:06:28.346351Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12, 12))\n",
    "ax = fig.add_subplot(1, 1, 1, projection='3d')\n",
    "\n",
    "ax.scatter3D(\n",
    "    image_rgb[:, :, 0],\n",
    "    image_rgb[:, :, 1],\n",
    "    image_rgb[:, :, 2],\n",
    "    c=image_rgb.reshape(-1, 3),\n",
    "    alpha=0.5,\n",
    ")\n",
    "ax.set_xlabel('red')\n",
    "ax.set_ylabel('green')\n",
    "ax.set_zlabel('blue')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:53.403058Z",
     "start_time": "2018-11-27T15:06:53.068733Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.utils import resample\n",
    "\n",
    "# perform k-means on a small random sample of pixels\n",
    "# we could use all pixels here but it would take much too long.\n",
    "# we reshape the image into a 2D array for resampling and k-means fitting\n",
    "flattened_image = image_hsv.reshape(-1, 3)\n",
    "sample = resample(flattened_image, n_samples=10000, replace=False)\n",
    "\n",
    "\n",
    "\n",
    "# get the desired number of cluster centers\n",
    "kmeans = KMeans(n_clusters = 50).fit(sample)\n",
    "centroids = kmeans.cluster_centers_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:53.840907Z",
     "start_time": "2018-11-27T15:06:53.405643Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(12, 4))\n",
    "\n",
    "centroids_rgb = color.hsv2rgb(centroids.reshape(1, -1, 3)).reshape(-1, 3)\n",
    "\n",
    "ax.scatter(\n",
    "    np.linspace(0, 1, len(centroids_rgb)),\n",
    "    np.ones(len(centroids_rgb)),\n",
    "    c=centroids_rgb,\n",
    "    s=10000.\n",
    ")\n",
    "ax.set_axis_off()\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:54.792692Z",
     "start_time": "2018-11-27T15:06:53.843275Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# associate each pixel with the number of a cluster (i.e. the nearest centroid/color)\n",
    "labels = kmeans.predict(flattened_image)\n",
    "\n",
    "# get the actual value for each cluster centroid and reshape into a 3D image\n",
    "reconstructed  = centroids[labels].reshape(width, length, channels)\n",
    "\n",
    "# convert to RGB and plot again.\n",
    "plt.imshow(color.hsv2rgb(reconstructed))\n",
    "plt.figure()\n",
    "plt.imshow(image_rgb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sorting colors using PCA\n",
    "\n",
    "When plotting the image centroids there are various ways to arrange the colors. There is no natural order \n",
    "to colors. Just as there is no natural order to most things with more than 1 dimension. \n",
    "\n",
    "We can try a few different things and try to get a nice output.\n",
    "\n",
    "1. Order by hue, Value, Saturation in HSV colorspace\n",
    "2. Order by red, green, blue in RGB colorspace\n",
    "3. Order by *length* of vector in e.g. eucledian norm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:55.375841Z",
     "start_time": "2018-11-27T15:06:54.795259Z"
    }
   },
   "outputs": [],
   "source": [
    "flattened_image = image_hsv.reshape(-1, 3)\n",
    "sample = resample(flattened_image, n_samples=10000, replace=False)\n",
    "\n",
    "# get the desired number of cluster centers\n",
    "kmeans = KMeans(n_clusters = 100).fit(sample)\n",
    "\n",
    "hsv = kmeans.cluster_centers_\n",
    "rgb = color.hsv2rgb(hsv.reshape(1, -1, 3)).reshape(-1, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:56.425540Z",
     "start_time": "2018-11-27T15:06:55.378062Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "fig, axs = plt.subplots(6, 1, figsize=(12, 12))\n",
    "\n",
    "\n",
    "print(len(hsv), len(rgb))\n",
    "\n",
    "for i, label, ax in zip(range(3), 'rgb', axs[:3]):\n",
    "    \n",
    "    idx = np.argsort(rgb[:, i])\n",
    "    \n",
    "    ax.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "    ax.set_axis_off()\n",
    "    ax.set_title(f'Order by {label}')\n",
    "\n",
    "    \n",
    "\n",
    "for i, label, ax in zip(range(3), 'hsv', axs[3:]):\n",
    "    \n",
    "    idx = np.argsort(hsv[:, i])\n",
    "    \n",
    "    ax.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "    ax.set_axis_off()\n",
    "    ax.set_title(f'Order by {label}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:57.328236Z",
     "start_time": "2018-11-27T15:06:56.937404Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 6))\n",
    "\n",
    "\n",
    "norm = np.linalg.norm(rgb, axis=1)\n",
    "idx = np.argsort(norm)\n",
    "\n",
    "ax1.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax1.axis('off')\n",
    "ax1.set_title('Order by norm RGB')\n",
    "\n",
    "norm = np.linalg.norm(centroids, axis=1)\n",
    "idx = np.argsort(norm)\n",
    "ax2.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax2.axis('off')\n",
    "ax2.set_title('Order by norm HSV')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need some other way to reduce the information down to 1D so we can sort it. \n",
    "Previously we used the PCA to do just that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:57.333572Z",
     "start_time": "2018-11-27T15:06:57.330440Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:57.768544Z",
     "start_time": "2018-11-27T15:06:57.335763Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2,) = plt.subplots(2, 1, figsize=(12, 6))\n",
    "\n",
    "c = PCA(n_components=1).fit_transform(rgb).ravel()\n",
    "idx = np.argsort(c)\n",
    "\n",
    "ax1.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax1.axis('off')\n",
    "ax1.set_title('Order by PCA transformed RGB.')\n",
    "\n",
    "\n",
    "c = PCA(n_components=1).fit_transform(hsv).ravel()\n",
    "idx = np.argsort(c)\n",
    "\n",
    "ax2.scatter(np.linspace(0, 1, len(idx)), np.ones_like(idx),  c=rgb[idx], s=10000)\n",
    "ax2.axis('off')\n",
    "ax2.set_title('Order by PCA transformed RGB')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density Based Clustering\n",
    "\n",
    "K-Means only takes distances between neighbours into acount. This works well for convex distributions that are approximately round in shape. Below you'll find another example where k-Means fails miserably. The moons are neither round nor convex yet easily seperated into clusters by eye. A simple transformation will not help in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:06:58.165627Z",
     "start_time": "2018-11-27T15:06:57.771441Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_moons\n",
    "\n",
    "np.random.seed(42)\n",
    "X, y = make_moons(n_samples=400, noise=0.1)\n",
    "\n",
    "prediction = KMeans(n_clusters=2, random_state=0).fit_predict(X)\n",
    "\n",
    "\n",
    "plt.scatter(X[:, 0], X[:, 1], c=prediction, s=200, lw=0, alpha=0.5, cmap=discrete_cmap)\n",
    "plt.scatter(X[:, 0], X[:, 1], c=y, cmap=discrete_cmap)\n",
    "plt.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DBSCAN\n",
    "\n",
    "Density-based spatial clustering of applications with noise (DBSCAN) is a clsutering method based on the density of different regions within the parameter space. It clsuters together points within regions of high density and labels points as outliers that do not lie within a dense region.\n",
    "\n",
    "1. Start by selecting *core points*. These are all points that have at least $m$ points in their neighbourhood region of radius $\\epsilon$.\n",
    "\n",
    "2. For each core point find connected points to build a cluster $C$. A point is connected to $C$ if its within the $\\epsilon$  neighbourhood of any point in $C$ and is also a core point.\n",
    "\n",
    "3. Assign each point which is not a core point to the nearest cluster with distance being at most $\\epsilon$ \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:00.159918Z",
     "start_time": "2018-11-27T15:06:58.167799Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_moons, make_checkerboard, make_circles\n",
    "from sklearn.cluster import DBSCAN\n",
    "\n",
    "discrete_cmap = ListedColormap([f'C{i}' for i in range(2)])\n",
    "discrete_cmap.set_under('k')\n",
    "\n",
    "f, [top, center, bottom] = plt.subplots(3, 4, sharex=True, sharey=True, figsize=(20, 20))\n",
    "\n",
    "X_moon, y_moon = make_moons(noise=0.05, random_state=0)\n",
    "X_circle, y_circle = make_circles(noise=0.05, factor=0.4, random_state=0, n_samples=200)\n",
    "X_blobs, y_blobs = make_blobs(centers=2, center_box=(-0.5, 0.5), cluster_std=0.4, random_state=0)\n",
    "X_long, y_long = make_blobs(centers=2, center_box=(-2.1, 2.1), cluster_std=0.1, random_state=0)\n",
    "\n",
    "data = [(X_moon, y_moon), (X_long, y_long), (X_circle, y_circle), (X_blobs, y_blobs)]\n",
    "\n",
    "for ax, (X, y) in zip(top, data):\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=y, cmap=discrete_cmap, vmin=-0.5, vmax=1.5)\n",
    "    ax.axis('off')\n",
    "    \n",
    "for ax, (X, y) in zip(center, data):\n",
    "    prediction = KMeans(n_clusters=2, random_state=0).fit_predict(X)\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=prediction, cmap=discrete_cmap, vmin=-0.5, vmax=1.5)\n",
    "    ax.axis('off')\n",
    "    \n",
    "for ax, (X, y) in zip(bottom, data):\n",
    "    prediction = DBSCAN(eps=0.339).fit_predict(X)\n",
    "    ax.scatter(X[:, 0], X[:, 1], c=prediction, cmap=discrete_cmap, vmin=-0.5, vmax=1.5)\n",
    "    ax.axis('off')\n",
    "    \n",
    "top[1].set_title('True Clusters')\n",
    "center[1].set_title('K - Mean Clusters')\n",
    "bottom[1].set_title('DBSCAN Clusters')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluating Cluster Performance\n",
    "\n",
    "In the case of unsupervised learning there is no ground truth to which the cluster structure can be compared. \n",
    "\n",
    "Some heuristic has to be applied measure how well a clustering performed. \n",
    "\n",
    "#### Silhuette Coefficent \n",
    "\n",
    "This evaluation ciriterion assumes a clsutering is 'good' if the clusters are dense instead of sparse. \n",
    "\n",
    "Define $a$ as the distance between a single point $x_0$ and all other points in the cluster $S_p$.\n",
    "\n",
    "$$\n",
    "a(x_0) = \\sum_{x_i \\in S_0} \\left\\|x_i - x_0 \\right\\|\n",
    "$$\n",
    "\n",
    "and $b$ as the distance between $x_0$ and all the points in the *nearest cluster* $S_p^\\prime$\n",
    "\n",
    "$$\n",
    "b(x_0) = \\sum_{x_i \\in S_0^\\prime} \\left\\|x_i - x_0 \\right\\|\n",
    "$$\n",
    "\n",
    "The Silhuette Coefficent is then defined as \n",
    "\n",
    "$$\n",
    "s = \\frac{b - a}{\\text{max}(a, b)}\n",
    "$$\n",
    "\n",
    "\n",
    "The coefficent takes a value close to +1 for dense clustering and -1 for sparse clusters. \n",
    "\n",
    "Unfortunately it doesn't work very reliable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:00.582110Z",
     "start_time": "2018-11-27T15:07:00.162684Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import silhouette_score\n",
    "\n",
    "np.random.seed(20)\n",
    "X, y = make_moons(n_samples=360, noise=0.1, random_state=172)\n",
    "\n",
    "km = KMeans(n_clusters=2)\n",
    "prediction_kmeans = km.fit_predict(X)\n",
    "score_kmeans = silhouette_score(X, km.labels_) \n",
    "\n",
    "db = DBSCAN(eps=0.14)\n",
    "prediction_db = db.fit_predict(X)\n",
    "score_db = silhouette_score(X, db.labels_) \n",
    "\n",
    "fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(10, 5))\n",
    "ax1.set_title('k-Means clustering score: {:0.3f}'.format(score_kmeans))\n",
    "ax1.scatter(X[:, 0], X[:, 1], c=prediction_kmeans, cmap=discrete_cmap, vmin=-0.5, vmax=1.5)\n",
    "ax1.axis('off')\n",
    "\n",
    "ax2.set_title('DBSCAN clustering: {:0.3f}'.format(score_db))\n",
    "ax2.scatter(X[:, 0], X[:, 1], c=prediction_db, cmap=discrete_cmap, vmin=-0.5, vmax=1.5)\n",
    "ax2.axis('off')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Complete Example\n",
    "\n",
    "Below we load a dataset containing data observed by the FACT telescope.\n",
    "\n",
    "<img width=\"45%\" src=\"http://www.miguelclaro.com/wp/wp-content/uploads/2013/10/FACTMilkyWayVertical-4650-net.jpg\" />   \n",
    "\n",
    "We will perform the typical steps to build and evaluate a classifier.\n",
    "\n",
    "0. Understand where your data comes from\n",
    "\n",
    "1. Preprocessing\n",
    "    \n",
    "    * Drop Constant Values,\n",
    "    * Handle Missing Data \n",
    "    * Feature Generation\n",
    "\n",
    "2. Splitting\n",
    "    \n",
    "    * Split your data into training and evaluation sets\n",
    "    \n",
    "3. Training \n",
    "    \n",
    "    * Train your classifier of choice.\n",
    "    \n",
    "4. Evaluation\n",
    "    \n",
    "    * Evaluate the performance on the test data set.\n",
    "    * If not good enough, go back to step 1 \n",
    "    \n",
    "5. Physics\n",
    "    \n",
    "    * Check whether your data support your hypothesis\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Get to know your data\n",
    "\n",
    "Cherenkov telescopes record short flashes of light produced by very high energy cosmic rays and photons hitting earths atmosphere.\n",
    "\n",
    "![](https://www.cta-observatory.org/wp-content/uploads/2016/05/cta47.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:00.589316Z",
     "start_time": "2018-11-27T15:07:00.584438Z"
    }
   },
   "outputs": [],
   "source": [
    "%%HTML\n",
    "\n",
    "<!-- https://nextcloud.e5.physik.tu-dortmund.de/index.php/s/e7yb2mifGDeyDBN/download -->\n",
    "<video width=\"100 %\" controls>\n",
    "  <source src=\"./resources/fact_events.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use machine learning for two tasks in this example. \n",
    "\n",
    " * Train a classifier to distinguish events induced by gamma rays form events induced by cosmic rays\n",
    " * Train a regressor to estimate the energy of the incoming primary particle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Preprocess data\n",
    "\n",
    "Load data and remove unwanted columns and store the true labels separately.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:01.600004Z",
     "start_time": "2018-11-27T15:07:00.592824Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import fact.io as fio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.004702Z",
     "start_time": "2018-11-27T15:07:01.603010Z"
    }
   },
   "outputs": [],
   "source": [
    "gammas = fio.read_h5py('./resources/sample_diffuse_gammas.h5', key='events')\n",
    "gammas.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now delete all simulated values which can not be observed during measurement in the physical world. We know which columns to remove because they have a special prefix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.055025Z",
     "start_time": "2018-11-27T15:07:02.007914Z"
    }
   },
   "outputs": [],
   "source": [
    "forbidden_columns = 'ceres_|mc_|corsika_|run_|source_position_|pointing_|aux_|event_num|incident_angle'\n",
    "\n",
    "gammas = gammas.filter(regex=f'^(?!{forbidden_columns}).*$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the data types of the columns. We can select non-numeric types and encode them. But in this case we might as well drop them as the attribute is not important."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.075296Z",
     "start_time": "2018-11-27T15:07:02.057255Z"
    }
   },
   "outputs": [],
   "source": [
    "c = gammas.select_dtypes(exclude=['number']).columns\n",
    "print(c)\n",
    "\n",
    "gammas = gammas.drop(c, axis='columns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can spot the columns with constant values by looking at the standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.529060Z",
     "start_time": "2018-11-27T15:07:02.078743Z"
    }
   },
   "outputs": [],
   "source": [
    "desc = gammas.describe()\n",
    "desc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.547778Z",
     "start_time": "2018-11-27T15:07:02.532415Z"
    }
   },
   "outputs": [],
   "source": [
    "c = desc.columns[desc.loc['std'] == 0]\n",
    "print(c)\n",
    "gammas = gammas.drop(c, axis='columns')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check for missing data. (Just delete it in this case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.594941Z",
     "start_time": "2018-11-27T15:07:02.551401Z"
    }
   },
   "outputs": [],
   "source": [
    "gammas = gammas.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we only loaded simulated gamma-ray showers. Now we do the same for the cosmic ray events. We create a method to perform all preprocessing in one step. We need this several times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:02.608429Z",
     "start_time": "2018-11-27T15:07:02.597366Z"
    }
   },
   "outputs": [],
   "source": [
    "def preprocess(df):\n",
    "    df = df.filter(regex=f'^(?!{forbidden_columns}).*$')\n",
    "    c = df.select_dtypes(exclude=['number']).columns\n",
    "    df = df.drop(c, axis='columns')\n",
    "    desc = df.describe()\n",
    "    c = desc.columns[desc.loc['std'] == 0]\n",
    "    df = df.drop(c, axis='columns')\n",
    "    df = df.dropna()\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:03.467654Z",
     "start_time": "2018-11-27T15:07:02.611649Z"
    }
   },
   "outputs": [],
   "source": [
    "protons = fio.read_h5py('./resources/sample_proton.h5', key='events')\n",
    "protons = preprocess(protons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can perform feature generation. We use our expert knowledge or intuition to create a new feature by combining existing columns into a new variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:03.481076Z",
     "start_time": "2018-11-27T15:07:03.471131Z"
    }
   },
   "outputs": [],
   "source": [
    "def feature_generation(df):\n",
    "    df['awesome_feature'] = (df.length/df.width) * df.size\n",
    "    return df\n",
    "\n",
    "gammas = feature_generation(gammas)\n",
    "protons = feature_generation(protons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A quick look at the data so far"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:03.829801Z",
     "start_time": "2018-11-27T15:07:03.484278Z"
    }
   },
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 20, 100)\n",
    "plt.hist(gammas.width, bins=bins, histtype='step', lw=4)\n",
    "plt.hist(protons.width, bins=bins, histtype='step', lw=4)\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we combine the two datasets into one big matrix and build a label vector $y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:03.895879Z",
     "start_time": "2018-11-27T15:07:03.832960Z"
    }
   },
   "outputs": [],
   "source": [
    "y = np.concatenate([np.ones(len(gammas)), np.zeros(len(protons))])\n",
    "X = pd.concat([gammas, protons])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Split Data\n",
    "\n",
    "Now we can split the data into test and training sets. Scikit-Learn provides some neat methods to do just that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:07:03.983514Z",
     "start_time": "2018-11-27T15:07:03.898708Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X_test, X_train, y_test, y_train = train_test_split(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Train the classifier\n",
    "\n",
    "Now we can train any classifier we want on the prepared data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:13:04.181685Z",
     "start_time": "2018-11-27T15:13:02.875532Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "rf = DecisionTreeClassifier(max_depth=10, criterion='entropy')\n",
    "rf.fit(X_train, y_train)\n",
    "\n",
    "y_prediction = rf.predict(X_test)\n",
    "y_prediction_proba = rf.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Evaluation \n",
    "\n",
    "Check accuracy of the models and other metrics "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:13:04.199205Z",
     "start_time": "2018-11-27T15:13:04.183884Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score, roc_curve, roc_auc_score\n",
    "\n",
    "acc = accuracy_score(y_test, y_prediction)\n",
    "auc = roc_auc_score(y_test, y_prediction_proba[:, 1])\n",
    "fpr, tpr, thresholds = roc_curve(y_test, y_prediction_proba[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:13:46.845086Z",
     "start_time": "2018-11-27T15:13:46.501506Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.scatter(fpr, tpr, c=thresholds)\n",
    "plt.plot(fpr, tpr, '--', color='gray', alpha=0.5)\n",
    "plt.text(0.5, 0.5, f'AuC ROC: {auc:0.03f} \\nAccuracy: {acc:0.03f}')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform steps 3, 4, and 5 in one step using cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:08:28.709336Z",
     "start_time": "2018-11-27T15:08:03.738680Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_validate\n",
    "\n",
    "rf = DecisionTreeClassifier(max_depth=12, criterion='entropy')\n",
    "\n",
    "scoring = {'acc': 'accuracy',\n",
    "           'auc': 'roc_auc',\n",
    "           'recall': 'recall'}\n",
    "\n",
    "results = cross_validate(rf, X, y, cv=5, scoring=scoring, return_train_score=True)\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:12:14.425262Z",
     "start_time": "2018-11-27T15:12:14.404548Z"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.display import Latex, display\n",
    "auc = results['test_auc']\n",
    "recall = results['test_recall']\n",
    "acc = results['test_acc']\n",
    "\n",
    "s = f'Area under RoC curve:  ${auc.mean():0.04f} \\\\pm {auc.std():0.04f}$'\n",
    "display(Latex(s))\n",
    "s = f' Accuracy:  ${acc.mean():0.04f} \\\\pm {acc.std():0.04f}$'\n",
    "display(Latex(s))\n",
    "s = f' Recall:  ${recall.mean():0.04f} \\\\pm {recall.std():0.04f}$'\n",
    "display(Latex(s))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Physics\n",
    "\n",
    "Now we could test our model and our hypothesis on real observed data. This part of the analysis is the most time \n",
    "consuming in general. It also requires more data than than this notebook can handle. \n",
    "After careful analysis one can produce an image of the gamma-ray sky\n",
    "\n",
    "<img width=\"60%\" src=\"https://www.mpi-hd.mpg.de/hfm/HESS/hgps/figures/HESS_J1813m126.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Improving Classification\n",
    "\n",
    "\n",
    "#### Boosting and AdaBoost\n",
    "\n",
    "Similar to the idea of combining many classifiers through bagging (like we did for the RandomForests) we now \n",
    "train many estimators in a sequential manner. In each iteration the data gets modified slightly using weights $w$\n",
    "for each sample in the training data. In the first iteration the weights are all set to $w=1$\n",
    "\n",
    "In each successive iteration the weights are updated. The samples that were incorrectly classified in the previous \n",
    "iteration get a higher weight. The weights for correctly classified samples get decreases. \n",
    "In other words: We increase the influence/importance of samples that are difficult to classify.\n",
    "\n",
    "Predictions are performed by taking a weighted average of the single predictors.\n",
    "\n",
    "The popular AdaBoost algorithms takes this a step further by optimizing the weight of each separate classifier \n",
    "in the ensemble.\n",
    "The AdaBoost ensemble combines many learners in an iterative way. The learner at iteration $m$ is:\n",
    "\n",
    "$$\n",
    " F_{m}(x)=F_{m-1}(x)+\\gamma _{m}h_{m}(x)\n",
    "$$\n",
    "\n",
    "The choice of $F_0$ is problem specific.\n",
    "\n",
    "Each weak learner produces a prediction $h(x_{m})$ for each sample in the training set. At each iteration $m$ a \n",
    "weak learner is fitted and assigned a coefficient $\\gamma_{m}$ which is found by minimizing:\n",
    "\n",
    "$$\n",
    "\\gamma_m = {\\underset {\\gamma }{\\arg \\min }} \\sum_{i}^{N}E\\bigl(F_{m-1}(x_{i})+\\gamma h(x_{i})\\bigr)\n",
    "$$\n",
    "\n",
    "where $E(F)$ is some error function and $x_i$ is the reweighted data sample.\n",
    "\n",
    "In general this method can work with any classifying method. Traditionally it is being used with very small \n",
    "decision trees. \n",
    "The weights get used to select the split points during the minimization of the loss function in each node\n",
    "\n",
    "$$\n",
    " \\underset{(X, s) \\in \\, \\mathbf{X} \\times {S}}{\\arg \\max} IG(X,Y) =   \\underset{(X, s) \\in \\, \\mathbf{X} \\times {S}}{\\arg \\max} ( H(Y) - H(Y |\\, X) ).\n",
    "$$\n",
    "\n",
    "Below we try AdaBoost on the FACT data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:49:09.247110Z",
     "start_time": "2018-11-27T15:48:58.872812Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import AdaBoostClassifier\n",
    "\n",
    "ada = AdaBoostClassifier(n_estimators=100, learning_rate=0.5, algorithm='SAMME')\n",
    "ada.fit(X_train, y_train)\n",
    "\n",
    "y_prediction = ada.predict(X_test)\n",
    "y_prediction_proba = ada.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:49:10.072785Z",
     "start_time": "2018-11-27T15:49:09.249195Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(range(100), list(ada.staged_score(X_test, y_test)), '.')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.xlabel('Iteration')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T15:49:10.727863Z",
     "start_time": "2018-11-27T15:49:10.075262Z"
    }
   },
   "outputs": [],
   "source": [
    "acc = accuracy_score(y_test, y_prediction)\n",
    "auc = roc_auc_score(y_test, y_prediction_proba[:, 1])\n",
    "fpr, tpr, thresholds = roc_curve(y_test, y_prediction_proba[:, 1])\n",
    "\n",
    "plt.scatter(fpr, tpr, c=thresholds)\n",
    "plt.plot(fpr, tpr, '--', color='gray', alpha=0.5)\n",
    "plt.text(0.5, 0.5, f'AuC ROC: {auc:0.03f} \\nAccuracy: {acc:0.03f}')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gradient Boosting \n",
    "\n",
    "Very similar to AdaBoost. Only this time we change the target label we train the classifiers for.\n",
    "\n",
    "Formulate the general problem as follows (See Wikipedia):\n",
    "\n",
    "Start with a constant function $F_{0}(x)$ and some differentiable loss function $L$ and incrementally expands it in a greedy fashion:\n",
    "\n",
    "$$\n",
    "F_{0}(x)={\\underset {\\gamma }{\\arg \\min }}{\\sum _{i=1}^{n}{L(y_{i},\\gamma )}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "F_{m}(x)=F_{m-1}(x)+{\\underset {h_{m}\\in {\\mathcal {H}}}{\\operatorname {arg\\,min} }}\\left[{\\sum _{i=1}^{n}{L(y_{i},F_{m-1}(x_{i})+h_{m}(x_{i}))}}\\right]\n",
    "$$\n",
    "\n",
    "Finding the best $ h_{m}\\in {\\mathcal {H}}$ is computationally speaking impossible.\n",
    "If we could find the perfect $h$ however, we know that \n",
    "\n",
    "$$\n",
    "F_{m+1}(x_i)=F_{m}(x_i)+h(x_i)=y_i\n",
    "$$\n",
    "\n",
    "or, equivalently, \n",
    "\n",
    "$$\n",
    "   h(x_i)= y_i - F_{m}(x_i)\n",
    "$$\n",
    "\n",
    "Note that for the mean squared error loss $\\frac{1}{2}(y_i - F(x_i))^2$ this is equivalent to the negative \n",
    "gradient with respect to $F_i$.\n",
    "\n",
    "For a general loss function we fit $h_{m}(x)$ to the residuals, or negative gradients \n",
    "$$\n",
    " r_{i, m}=-\\left[{\\frac {\\partial L(y_{i},F(x_{i}))}{\\partial F(x_{i})}}\\right]_{F(x)=F_{m-1}(x)}\\quad {\\mbox{for }}i=1,\\ldots ,n.\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "Below we try it on FACT data again.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T16:56:38.715276Z",
     "start_time": "2018-11-27T16:56:29.657159Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "\n",
    "grb = GradientBoostingClassifier(verbose=True, n_estimators=100)\n",
    "grb.fit(X_train, y_train)\n",
    "\n",
    "y_prediction = grb.predict(X_test)\n",
    "y_prediction_proba = grb.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T16:56:39.343691Z",
     "start_time": "2018-11-27T16:56:38.718176Z"
    }
   },
   "outputs": [],
   "source": [
    "l = [accuracy_score(y_pred, y_test) for y_pred in grb.staged_predict(X_test)]\n",
    "\n",
    "plt.plot(range(len(l)), l, '.')\n",
    "plt.ylabel('Accuracy')\n",
    "plt.xlabel('Iteration')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T16:56:40.215880Z",
     "start_time": "2018-11-27T16:56:39.346460Z"
    }
   },
   "outputs": [],
   "source": [
    "acc = accuracy_score(y_test, y_prediction)\n",
    "auc = roc_auc_score(y_test, y_prediction_proba[:, 1])\n",
    "fpr, tpr, thresholds = roc_curve(y_test, y_prediction_proba[:, 1])\n",
    "\n",
    "plt.scatter(fpr, tpr, c=thresholds)\n",
    "plt.plot(fpr, tpr, '--', color='gray', alpha=0.5)\n",
    "plt.text(0.5, 0.5, f'AuC ROC: {auc:0.03f} \\nAccuracy: {acc:0.03f}')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More on gradient descent algorithms can be found in the Neural Network notebook.\n",
    "\n",
    "Let's now test our all time favorite classifier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T16:58:57.873659Z",
     "start_time": "2018-11-27T16:58:45.510336Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "rf = RandomForestClassifier(n_estimators=150,  max_depth=18, criterion='entropy')\n",
    "rf.fit(X_train, y_train)\n",
    "\n",
    "y_prediction = rf.predict(X_test)\n",
    "y_prediction_proba = rf.predict_proba(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-27T16:58:58.442111Z",
     "start_time": "2018-11-27T16:58:57.875736Z"
    }
   },
   "outputs": [],
   "source": [
    "acc = accuracy_score(y_test, y_prediction)\n",
    "auc = roc_auc_score(y_test, y_prediction_proba[:, 1])\n",
    "fpr, tpr, thresholds = roc_curve(y_test, y_prediction_proba[:, 1])\n",
    "\n",
    "plt.scatter(fpr, tpr, c=thresholds)\n",
    "plt.plot(fpr, tpr, '--', color='gray', alpha=0.5)\n",
    "plt.text(0.5, 0.5, f'AuC ROC: {auc:0.03f} \\nAccuracy: {acc:0.03f}')\n",
    "None"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}