{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Tut5 - K_means_from_scratch.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EN6uJPFmZjnD"
      },
      "source": [
        "#K-means clustering"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "40Z7SQXeDujI"
      },
      "source": [
        "Kmeans clustering is one of the most popular clustering algorithms and usually the first thing practitioners apply when solving clustering tasks to get an idea of the structure of the dataset. Given a set of observations (x1, x2, ..., xn), where each observation is a d-dimensional real vector, k-means clustering aims to partition the n observations into k (≤ n) sets.\n",
        "\n",
        "Why do we cluster?\n",
        "* Exploratory data analysis technique used to get an intuition about the structure of the data.\n",
        "* Summarize data - Represent a large continuous vector with a cluster number.\n",
        "\n",
        "Algorithm:\n",
        "* Assume K clusters\n",
        "* Assume points that are close to other points belong to the same group. Far points are outside of the group\n",
        "* Randomly initialize cluster centers\n",
        "* Calculage distance of points to centers\n",
        "* **Assign** each point to closest cluster center\n",
        "* **Update**, cluster centers: center of gravity of points in each class.\n",
        "* Repeat!\n",
        "\n",
        "**Assignment step:** Assign each observation to the cluster with the nearest mean: that with the least squared Euclidean distance. \n",
        "\n",
        "![alt text](https://wikimedia.org/api/rest_v1/media/math/render/svg/602c759bad7686552299b5c1ecc3ff9188ced9b4)\n",
        "\n",
        "**Update step:** Recalculate means (centroids) for observations assigned to each cluster.\n",
        "\n",
        "![alt text](https://wikimedia.org/api/rest_v1/media/math/render/svg/964165b4792e24484c3ccde9b3e00c478c96fbd6)\n",
        "\n",
        "Output: \n",
        "  * centroids for each cluster\n",
        "  * cluster id for each point\n",
        "\n",
        "\n",
        "Design Choices:\n",
        "* Initialization\n",
        "  * Randomly select 𝐾 points as initial cluster centers\n",
        "  * Number of hyperparrameter 𝐾 is important\n",
        "\n",
        "Challenges:\n",
        "* Optimization\n",
        "  * Converges to a local minimum\n",
        "  * May want to perform multiple restarts (re-initialize and try again)\n",
        "\n",
        "![alt text](https://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/K-means_convergence.gif/220px-K-means_convergence.gif)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "s53FeZKzPUmi"
      },
      "source": [
        "from numpy.linalg import norm\n",
        "\n",
        "#Objecct Oriented Implementation of K-means\n",
        "class Kmeans:\n",
        "\n",
        "  #Set some parameters to initialize \n",
        "  def __init__(self, n_clusters, max_iter=100, tolerance = 1e-3,debug=False,random_state=123):\n",
        "    #nSpecify umber of centroids/clusters K\n",
        "    self.n_clusters = n_clusters\n",
        "    #Maximum iterations to run the algorithm\n",
        "    self.max_iter = max_iter\n",
        "    #Initialize random number generator\n",
        "    self.random_state = random_state\n",
        "    # Tolerance value to indicate that no improvement is made\n",
        "    self.tol = tolerance\n",
        "    # enable/disable debug messages\n",
        "    self.debug = debug\n",
        "  \n",
        "  #Initialize the centroids\n",
        "  def initialize_centroids(self, X):\n",
        "    #Initialize centroids by first shuffling the dataset and then randomly \n",
        "    #selecting K data points for the centroids without replacement\n",
        "    #np.random.RandomState(self.random_state)\n",
        "    random_idx = np.random.permutation(X.shape[0])\n",
        "    #Randomly initialize K cluster centroids\n",
        "    centroids = np.array(X[random_idx[:self.n_clusters]])\n",
        "    return centroids\n",
        "\n",
        "  #recompute a centroid based on the points x_i that belong to the same cluster.\n",
        "  def compute_centroids(self, X, labels):\n",
        "    centroids = np.zeros((self.n_clusters, X.shape[1]))\n",
        "    for k in range(self.n_clusters):\n",
        "      # average of points in the same cluster\n",
        "      centroids[k, :] = np.mean(X[labels == k, :], axis=0)\n",
        "    return centroids\n",
        "\n",
        "  #find distance of each point x_i from all the centroids\n",
        "  def compute_distance(self, X, centroids):\n",
        "    distance = np.zeros((X.shape[0], self.n_clusters))\n",
        "    for k in range(self.n_clusters):\n",
        "      row_norm = norm(X - centroids[k, :], axis=1)\n",
        "      distance[:, k] = np.square(row_norm)\n",
        "    return distance\n",
        "\n",
        "  #select the centroid with the smaller distance for each point\n",
        "  def find_closest_cluster(self, distance):\n",
        "    return np.argmin(distance, axis=1)\n",
        "\n",
        "  def compute_sse(self, X, labels, centroids):\n",
        "    distance = np.zeros(X.shape[0])\n",
        "    for k in range(self.n_clusters):\n",
        "      # Average (mean) of points assigned to cluster k\n",
        "      distance[labels == k] = norm(X[labels == k] - centroids[k], axis=1)\n",
        "    return np.sum(np.square(distance))    \n",
        "\n",
        "  # Main steps of K-Means\n",
        "  def fit(self, X):\n",
        "    self.centroids = self.initialize_centroids(X)\n",
        "    if(self.debug):\n",
        "      print('Initial centroids:')\n",
        "      print(self.centroids)\n",
        "\n",
        "    # track lowest error\n",
        "    min_error = np.Inf  \n",
        "\n",
        "    for i in range(self.max_iter):\n",
        "      # save previous centroids to compare how they change\n",
        "      old_centroids = self.centroids\n",
        "\n",
        "      #Compute the sum of the squared distance between data points and \n",
        "      #all centroids.\n",
        "      #find_distance of points x_i \n",
        "      distance = self.compute_distance(X, old_centroids)\n",
        "      # Assign each data point to the closest cluster (centroid).\n",
        "      # find index of cluster centroid closest to x_i\n",
        "      self.labels = self.find_closest_cluster(distance)\n",
        "      # Compute the centroids for the clusters by taking the average of the \n",
        "      # all data points that belong to each cluster.\n",
        "      # average (mean) of points assigned to cluster k\n",
        "      self.centroids = self.compute_centroids(X, self.labels)\n",
        "\n",
        "      #Compute cost function (distortion)\n",
        "      #distortion, is defined as the sum of the squared distances between each observation vector and its dominating centroid.\n",
        "      self.error = (1/X.shape[0])*self.compute_sse(X, self.labels, self.centroids)\n",
        "\n",
        "      # keep only the best centroids\n",
        "      if(self.error < min_error):\n",
        "        min_error = self.error\n",
        "        best_centroids = self.centroids\n",
        "\n",
        "      if(self.debug):\n",
        "        print(i,\": \",self.error)\n",
        "\n",
        "      # Keep iterating until there is no change to the centroids. i.e assignment \n",
        "      # of data points to clusters isn’t changing.\n",
        "      # If no significant change to any of the clusters then terminate  \n",
        "      if self.compute_distance(self.centroids, old_centroids).all()< self.tol:\n",
        "        # Pick clustering solution that gave the lowest distortion\n",
        "        self.centroids = best_centroids\n",
        "        break\n",
        "\n",
        "  def predict(self, X):\n",
        "    # find disances of x_i from centroids\n",
        "    distance = self.compute_distance(X, self.centroids)\n",
        "    #return the index of the closest centroid\n",
        "    return self.find_closest_cluster(distance)\n",
        "  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DccWXA9F5aiB"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "from matplotlib import style\n",
        "import matplotlib.cm as cm\n",
        "\n",
        "style.use('ggplot')\n",
        "\n",
        "#1.\n",
        "X = np.array([[1, 2],\n",
        "              [1.5, 1.8],\n",
        "              [5, 8 ],\n",
        "              [8, 8],\n",
        "              [1, 0.6],\n",
        "              [9,11]])\n",
        "y=None\n",
        "\n",
        "from sklearn.datasets import make_blobs,make_moons, load_iris\n",
        "# generate 2d classification dataset\n",
        "\n",
        "#2.\n",
        "#X, y = make_blobs(n_samples=100, centers=2, n_features=2)\n",
        "\n",
        "#3.\n",
        "#X, y = make_moons(n_samples=100, noise=0.1)\n",
        "\n",
        "#4.\n",
        "#This data sets consists of 3 different types of irises’ (Setosa, Versicolour, and Virginica) petal and sepal length,\n",
        "# iris = load_iris()\n",
        "# X = iris.data[:, :2]  # we only take the first two features.\n",
        "# y = iris.target\n",
        "\n",
        "print('Number of samples: ',X.shape[0])\n",
        "print('Number features: ',X.shape[1])\n",
        "plt.scatter(X[:, 0], X[:, 1],c=y)\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0Fsyf9R95k5C"
      },
      "source": [
        "# Run local implementation of kmeans\n",
        "\n",
        "model = Kmeans(n_clusters=2, max_iter=100,debug=False)\n",
        "model.fit(X)\n",
        "\n",
        "# Plot the clustered data\n",
        "fig, ax = plt.subplots(figsize=(6, 6))\n",
        "labnames=np.array([i+1 for i in range(model.n_clusters)]).T\n",
        "plt.scatter(X[:, 0], X[:, 1],c=model.labels)\n",
        "plt.scatter(model.centroids[:, 0], model.centroids[:, 1], marker='*', s=300, c='r', label='centroid')\n",
        "plt.title('Visualization of clustered data', fontweight='bold')\n",
        "ax.set_aspect('equal')\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7A3gXBXSIvXT"
      },
      "source": [
        "##How to choose the number of clusters?\n",
        "\n",
        "Elbow method: Stop adding clusters when improvement is small.\n",
        "\n",
        "Sometimes it’s still hard to figure out a good number of clusters to use because the curve is monotonically decreasing and may not show any elbow or has an obvious point where the curve starts flattening out."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V_wN8tnzBGgT"
      },
      "source": [
        "#Finding the best number of clusters\n",
        "X, y = make_blobs(n_samples=100, centers=3, n_features=2)\n",
        "X, y = make_moons(n_samples=100, noise=0.1)\n",
        "\n",
        "err = []\n",
        "for k in range(1,21):\n",
        "  me=np.Inf\n",
        "  for j in range(10):\n",
        "    model = Kmeans(n_clusters=k)\n",
        "    model.fit(X)\n",
        "    if(model.error<me):\n",
        "      me = model.error\n",
        "  err.append(me)\n",
        "plt.plot(err)\n",
        "plt.xticks(range(20),[r for r in range(1,21)])\n",
        "plt.xlabel('K - Number of clusters')\n",
        "plt.ylabel('Error')\n",
        "print(np.argmin(err))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RJWr6KKQZqR1"
      },
      "source": [
        "#K-Means and Spectral Clustering"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yLe6LcxaEwFk"
      },
      "source": [
        "K-Means suffers as the geometric shapes of clusters deviates from spherical shapes. The idea is we transform to higher dimensional representation that make the data linearly separable (the same idea that we use in SVMs). Different kinds of algorithms work very well in such scenarios such as SpectralClustering.\n",
        "\n",
        "Spectral clustering is a technique for clustering based on transformation by dimensionality reduction. In practice Spectral Clustering is very useful when the structure of the individual clusters is highly non-convex or more generally when a measure of the center and spread of the cluster is not a suitable description of the complete cluster. For instance when clusters are nested circles on the 2D plane.\n",
        "\n",
        "![alt text](https://image.slidesharecdn.com/spectralclustering-150305200247-conversion-gate01/95/spectral-clustering-12-638.jpg?cb=1425585864)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UlzURoxGHNnP"
      },
      "source": [
        "from sklearn.cluster import SpectralClustering\n",
        "model = SpectralClustering(n_clusters=2, affinity='nearest_neighbors',\n",
        "                           assign_labels='kmeans')\n",
        "labels = model.fit_predict(X)\n",
        "plt.scatter(X[:, 0], X[:, 1], c=labels,\n",
        "            s=50, cmap='viridis');"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BY5ti-YNZvWc"
      },
      "source": [
        "#K-Means for Image compression"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a7urJlPbJWbp"
      },
      "source": [
        "In computer graphics, color quantization is the task of reducing the color palette of an image to a fixed number of colors k.\n",
        "A color is represented as a three dimensional vector [Red, Green, Blue] each ranging from 0...255. So we have a total of 16581375 (255^3) colors. \n",
        "Our goal is to reduce the number of colors to k and represent (compress) the photo using those k colors only. \n",
        "we’ll use kmeans algorithm on the image and treat every pixel as a data point. That means reshape the image from height x width x channels to (height * width) x channel, i,e we would have 396 x 396 = 156,816 data points in 3-dimensional space which are the intensity of RGB. If we add also the x,y location we have (height * width) x (channel+2), i,e we would have 396 x 396 = 156,816 data points in 5-dimensional space \n",
        "\n",
        "![alt text](https://drive.google.com/uc?id=1Mo-l-xNpnk6z1DM9KUBjnelhTX_VZp4v)\n",
        "\n",
        "We do this with K-means by finding similar pixels and using the color corresponding to their centroid to represent them. \n",
        "\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0RK8pLgxMfw6"
      },
      "source": [
        "from skimage import data\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from PIL import Image\n",
        "from sklearn import preprocessing\n",
        "import cv2\n",
        "data.astronaut\n",
        "image = data.astronaut()\n",
        "image = cv2.resize(image,(100,100))\n",
        "plt.imshow(image)\n",
        "\n",
        "print(image.shape)\n",
        "imageH,imageW = image.shape[0:2]\n",
        "\n",
        "#\tInitialise data vector with attribute r,g,b,x,y for each pixel\n",
        "dataVector = np.ndarray(shape=(imageW * imageH, 5), dtype=float)\n",
        "#\tInitialise vector that holds which cluster a pixel is currently in\n",
        "pixelClusterAppartenance = np.ndarray(shape=(imageW * imageH), dtype=int)\n",
        "\n",
        "#\tPopulate data vector with data from input image\n",
        "#\tdataVector has 5 fields: red, green, blue, x coord, y coord\n",
        "for y in range(0, imageH):\n",
        "      for x in range(0, imageW):\n",
        "      \txy = (x, y)\n",
        "      \trgb = image[y,x,:]\n",
        "      \tdataVector[x + y * imageW, 0] = rgb[0]\n",
        "      \tdataVector[x + y * imageW, 1] = rgb[1]\n",
        "      \tdataVector[x + y * imageW, 2] = rgb[2]\n",
        "      \tdataVector[x + y * imageW, 3] = x\n",
        "      \tdataVector[x + y * imageW, 4] = y\n",
        "\n",
        "#Features of different types may have different scales.\n",
        "#For example, pixel coordinates on a 100×100 image vs. RGB color values in the range [0,1].\n",
        "#Problem: Features with larger scales dominate clustering.\n",
        "#Solution: Scale the features.\n",
        "\n",
        "#\tStandarize the values of our features\n",
        "# RGB values and dimensions have different scales\n",
        "dataVector_scaled = preprocessing.normalize(dataVector)\n",
        "\n",
        "print(dataVector_scaled.shape)\n",
        "\n",
        "model = Kmeans(n_clusters=10,debug=False)\n",
        "model.fit(dataVector_scaled)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ro-QcjQitlYy"
      },
      "source": [
        "#Find the closest color for each datapoint\n",
        "import matplotlib.cm as cm\n",
        "colors = cm.rainbow(np.linspace(0, 1, len(model.centroids)))\n",
        "for ndv,dV in zip(dataVector_scaled,dataVector):\n",
        "  item = int(model.predict(ndv[None,:]))\n",
        "  dV[0] = int(round(model.centroids[item][0] * 255))\n",
        "  dV[1] = int(round(model.centroids[item][1] * 255))\n",
        "  dV[2] = int(round(model.centroids[item][2] * 255))\n",
        "\n",
        "# Replace all the pixels in the original image with the centroid color\n",
        "for y in range(imageH):\n",
        "\tfor x in range(imageW):\n",
        "\t \timage[y, x,:]=[int(dataVector[y * imageW + x][0]), \n",
        "\t \t\t\t\t\t\t\tint(dataVector[y * imageW + x][1]),\n",
        "\t \t\t\t\t\t\t\tint(dataVector[y * imageW + x][2])]\n",
        "\n",
        "#display the image\n",
        "plt.imshow(image)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kmhi1yoPPeBr"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}