L11.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "L11.ipynb",
      "provenance": []
    },
    "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.8"
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oJVOsD5TDbPK"
      },
      "source": [
        "In this notebook, we will write our own k-means model from scratch and use it to cluster handwritten numbers from the MNIST dataset\n",
        "\n",
        "\n",
        "In the cell below, we create an `assign_data` function, which takes the `data` and the `centers` for each cluster and makes an assignment of each datapoint in `data` to the closest of the `centers`, `centerids`.  We extract `n`, the number of datapoints, `d`, the dimensionality of the datapoints, and `k` the number of centers.\n",
        "\n",
        "Next, we need to compute the squared distance between each center and each data point.\n",
        "\n",
        "Reshaping the data to be 1 x `n` x `d`, and the centers to be `k` x 1 x `d` signals to numpy that when it subtracts these two arrays, it creates an array of shape `k` x `n` x `d`. That is, it computes all combinations of the `k` centers and the `n` datapoints for each of the `d` dimensions. We assign those differences to `res`.\n",
        "\n",
        "Squaring each of the differences, then summing along dimension 2 --- that’s the components of the vectors --- produces the sum of squared distances, which is the squared distance between the centers and the datapoints.  The resulting array is of shape `k` x `n`.\n",
        "\n",
        "`assign_data` also computes the `loss`, which the sum of the squared differences. We want to know which center has the *smallest* squared distance for each data point. `argmin` produces the index of an array with the smallest value along the given dimension. Here, we’re using dimension 0, which varies over the `k` centers. `centerids` is now an array with one integer for each datapoint that indicates which of the centers is closest.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QtZgh0QsJpid"
      },
      "source": [
        "def assign_data(data,centers):\n",
        "  # n is the number of data points\n",
        "  n = len(data)\n",
        "  # d is the dimensionality of the data points\n",
        "  d = len(data[0])\n",
        "  # k is the number of clusters\n",
        "  k = len(centers)\n",
        "  # first, subtract the set of centers from each data point\n",
        "  res = np.reshape(data,(1,n,d))-np.reshape(centers,(k,1,d))\n",
        "  # sum the squared differences\n",
        "  res2 = np.add.reduce(res**2,2)\n",
        "  # assign each data point to its closest center\n",
        "  centerids = np.apply_along_axis(np.argmin,0,res2)\n",
        "  # While we're here, make a note of the loss\n",
        "  loss = sum(np.apply_along_axis(np.min,0,res2))\n",
        "  return(centerids, loss)"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2x4tp5vkDqDh"
      },
      "source": [
        "Next we'll compute the mean of each of the `k` centers using the `data` and their `centerids` assignments. `compute_means` takes the data and the center ids and computes the centers by averaging all of the datapoints with the same id.  This will be used to update the `centers`.\n",
        "\n",
        "After extracting the number of datapoints and dimension, we initialize the array of center locations to a `k` x `d` array of all zeros.\n",
        "\n",
        "For each of the cluster id values from 0 to `k`, we do the following operations:\n",
        " - First, form a smaller array, `cols`, consisting of all the datapoints with the current center id.\n",
        " - To be robust, we make sure `cols` has a length greater than zero. That can happen if there’s a center that has been elbowed out of the running by the other centers being closer to all of the data points.\n",
        "  - If it equals zero, that means our center is out of the action and we should probably pick a different location for it. We simply choose one of the data points at random to be this new location.\n",
        "\n",
        "- We want to move the center to the `mean` of the closest points. Numpy’s `mean` method computes the average of an array along any given dimension. Here, we choose dimension 0, which corresponds to the different data points. `mean` produces a component-wise average of all the data points with cluster id equal to `i`.\n",
        "- After completing the loop, we return the newly computed `centers`.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vifay3EeJqi-"
      },
      "source": [
        "def compute_means(data, centerids, k):\n",
        "  # n is number of data points\n",
        "  n = len(data)\n",
        "  # d is dimensionality of the data points\n",
        "  d = len(data[0])\n",
        "\n",
        "  # Zero out the centers\n",
        "  centers = np.zeros(shape=(k,d))\n",
        "\n",
        "  # loop over the clusters\n",
        "  for i in range(k):\n",
        "    # Gather the data points assigned to cluster i\n",
        "    cols = np.array([data[j] for j in range(n) if centerids[j] == i])\n",
        "    # Average to get mean for that cluster\n",
        "    if len(cols) == 0:\n",
        "      centers[i] = data[random.randint(0,n-1)]\n",
        "    else:\n",
        "      centers[i] = cols.mean(0)\n",
        "  return(centers)"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5-qbkig5D29K"
      },
      "source": [
        "With these two functions, we can build a `kmeans` model, which takes in the `data` and number of clusters, `k` and iteratively builds `k` clusters and updates them relative to the `loss`.  \n",
        "\n",
        "We initialize the `k` centers by selecting random data points. We loop until the `loss` stops changing. If `oldloss` is different from the new `loss`, we use `assign_data` to assign each datapoint to its closest center. Then, we use `compute_means` to move the centers to the means of the points assigned to them. We repeat until the `loss` stops changing, returning the final `loss` and `centers`.  \n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ex_bJQGsJql8"
      },
      "source": [
        "def kmeans(data, k):\n",
        "  n = len(data)\n",
        "  d = len(data[0])\n",
        "  # grab the centers from random points\n",
        "  centers = data[[random.randint(0,n-1) for i in range(k)]]\n",
        "  oldloss = 0\n",
        "  loss = 1\n",
        "  while oldloss != loss:\n",
        "    oldloss = loss\n",
        "    centerids, loss = assign_data(data,centers)\n",
        "    centers = compute_means(data, centerids, k)\n",
        "  return(centers, loss)"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zxwnWjHWD7Cd"
      },
      "source": [
        "We will download the MNIST dataset and split the data into training data, `X_train` and `y_train` and test data, `X_test` and `y_test`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DZohJqLrJqoS",
        "outputId": "74394c41-d364-4ff9-f2ed-d1beeaa50fcf",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "source": [
        "from sklearn.datasets import fetch_openml\n",
        "# data = fetch_openml(name='mnist_784')\n",
        "data = fetch_openml(name='mnist_784', as_frame=False) # Thanks to SpencerGreene for this fix! (9/17/23)\n",
        "\n",
        "import numpy as np\n",
        "from sklearn.model_selection import train_test_split\n",
        "X_train, X_test, y_train, y_test = train_test_split(data.data, data.target, test_size=0.1)\n",
        "X_train, X_test, y_train, y_test = train_test_split(X_test, y_test, test_size=0.33)\n",
        "len(X_train)"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/dist-packages/sklearn/datasets/_openml.py:968: FutureWarning: The default value of `parser` will change from `'liac-arff'` to `'auto'` in 1.4. You can set `parser='auto'` to silence this warning. Therefore, an `ImportError` will be raised from 1.4 if the dataset is dense and pandas is not installed. Note that the pandas parser may return different data types. See the Notes Section in fetch_openml's API doc for details.\n",
            "  warn(\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "4690"
            ]
          },
          "metadata": {},
          "execution_count": 5
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NsiSIN2rELNf"
      },
      "source": [
        "Here, we run `kmeans` on our `X_train` data where `k=10`.  We run `kmeans` 9 times and record the `bestcenters` which have the `bestloss` among the recorded losses.  We then find the accuracy of these new centers on the test set."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "85T_3_oHKa31",
        "outputId": "858cfbc0-6237-43bc-8959-e8320bb42c83",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "source": [
        "from scipy import stats\n",
        "import math\n",
        "from functools import reduce\n",
        "import random\n",
        "\n",
        "#for nlabeled in range(20,len(X_train),10):\n",
        "nlabeled = 20\n",
        "if True:\n",
        "  print(nlabeled)\n",
        "  ans = []\n",
        "  k = 10 # 2 # 5 # 20\n",
        "  if True:\n",
        "    bestcenters, bestloss = kmeans(X_train, k)\n",
        "    for rep in range(9):\n",
        "      centers, loss = kmeans(X_train, k)\n",
        "      if loss < bestloss: bestcenters, bestloss = centers, loss\n",
        "    # How do we test the clustering that was discovered?\n",
        "    # Assign testing points to clusters\n",
        "    test_centerids, loss = assign_data(X_test, bestcenters)\n",
        "\n",
        "    ## Use the labeled examples to label the clusters\n",
        "    # First, assign the data points to clusters\n",
        "    train_centerids, loss = assign_data(X_train[:nlabeled], bestcenters)\n",
        "    # Now, look up the actual labels for each of the data points\n",
        "    labs = y_train[:nlabeled]\n",
        "\n",
        "    # Here's where we will put the label for each of the k clusters\n",
        "    clust_labs = np.repeat(labs[0],k)\n",
        "    # Run through each of the clusters, i\n",
        "    for i in range(k):\n",
        "      # Look up all the labels assigned to all the data points in cluster i, returning the most common one\n",
        "      mode = stats.mode(labs[train_centerids == i].astype(int)).mode # fix by SpencerGreene! (9/17/23)\n",
        "      if not np.isnan(mode): clust_labs[i] = mode # ditto\n",
        "\n",
        "    ans = ans + [(k,sum(clust_labs[test_centerids] == y_test)/len(y_test))]\n",
        "#    plt.plot(X_test[clust_labs[test_centerids] == 0,0],X_test[clust_labs[test_centerids] == 0,1],'o',color='r')\n",
        "#    plt.plot(X_test[clust_labs[test_centerids] == 1,0],X_test[clust_labs[test_centerids] == 1,1],'o',color='b')\n",
        "#    plt.show()\n",
        "\n",
        "#  print(ans)\n",
        "  print(reduce((lambda x, y: x if x[1] > y[1] else y), ans))\n",
        "  labids, loss = assign_data(X_test, X_train[:nlabeled])\n",
        "  print(nlabeled, sum(y_train[labids] == y_test)/len(y_test))"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "20\n",
            "(10, 0.3476190476190476)\n",
            "20 0.45151515151515154\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3iz6oNKK8Tux"
      },
      "source": [
        "We'll next print the images that best represent the centers of each of our clusters in K-means and the label for each of the clusters\n",
        "\n",
        "We will also calcuate the percent accuracy of the clusters"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "er8dfhHEJqrA",
        "outputId": "e45bea51-cbca-4b4e-d75e-833ca7fcd531",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 493
        }
      },
      "source": [
        "# !pip install keras=='2.8.0'\n",
        "from tensorflow import keras\n",
        "\n",
        "train_centerids, loss = assign_data(X_train, bestcenters)\n",
        "test_centerids, loss = assign_data(X_test, bestcenters)\n",
        "\n",
        "clust_labs = np.repeat(labs[0],k)\n",
        "for i in range(len(bestcenters)):\n",
        "  display(keras.utils.array_to_img(np.reshape(bestcenters[i],(28,28,1)), scale=False))\n",
        "  clust_labs[i] = y_train[train_centerids == i][0]\n",
        "  print(clust_labs[i])\n",
        "#  mode = stats.mode(y_train[train_centerids == i]).mode\n",
        "#  print(mode[0])\n",
        "#  if len(mode) > 0: clust_labs[i] = mode[0]\n",
        "\n",
        "sum(clust_labs[test_centerids] == y_test)/len(y_test)"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABp0lEQVR4nLWSW2/TQBCFZ/Zir6kT2ySpExqiVqRtoAr8/5/AC6AKCUhLQCglMcRN4pCLL+vhoTgxVcUb87b76ZydnTkA/yh88IoegMgMw5Z6nSSZvgcRuar7j6v59HYSbYkAxB89ISIazdOTDoUYMcAcCkiACMDdzkXXC2/C5SYjAABWmBKB2broHcXj0Y9ZrP+CRMCcZ712cv3+2zwlRNzZAgCgOnnVxavL4SplCIBUsgVx+LLvjS6/zFLBkIDKb7LKWf/p7ZvBzxgFx7sfFrYoW/0e/zQIQSmuEYH2ELlzfF6/+Ti1au7BZrJOdFlp1I476TDwGocOC8Vym5Ugf+S3zfC70X3i2KYb2bNtqVtZqbsw3dZqSqeqXlVS7JXMcDwrjqTS46TlSCQGuBu8NCqOyHMZTdeWeTBfx0m+Gx+CEEoy2wTWePHcCoJFXO6W0s3K96HCGj7/fPV1mQEgFSuLF0HryHVzzlbXr9+O1xntlk06DoeWeerxX+PBuw+Tlc6B9jHhht1sN6tiPhkFiyS/FzDGpTIYJptU6yJ8pYAh0N2R4D/Xb1SarmEEuTjfAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "6\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABjUlEQVR4nL3S23LaMBAG4H9XPgiMS1wKmXKTzrST93+f9KplpikBUscYHyTt9gKbQB+gupL0ze7srgT8/0XvO1IAIAZEoVdIZMgwi5CJyTsfAgBEA8U2yWaJISaty7pxFPSMFEV2XnxaFVYQyG9/vogEgiICmGP7Yf3wbYnj0SW5zUPXtBjTkkmL9eMX2j7v2unXlS62RkXHtJzYu5kvn37sw3JF0zxVER0Lgkkmptp835wi62gS4VwsIkBVJWJpn/cV0snMKvUO7yiibNjYjzxfL/LGO70gFKFvg/0sfZrd3887L8NsIgDq+/pQFMupEuVzC+dALEO1CF25MU0G9ZIY+GPjARpQyJVy+p0Z9Pzw2O0OrReVsRVxtau2hoFsEY77qlPVy+A1BE8VERlOUL31crZhCApPRGCaZfpWB71BAAoQ7N2k/9MKD6/MuFKT5UlTdWCi20gCKJmmqOpOxrsRiYhNHIWXw+711AW5QQUU4fSr5vK1bLxef7BzdGxTI6Hvgwj+RRCIVfRy/gsrmdQMnKnmXQAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "7\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABsUlEQVR4nLXSTY/TMBAG4BnbsZO4TdJvVpSyF0D7/3/J3hCrBbYLatmqaprWbRLbGQ5toy5I3JiTrUczmpFegH8UXr8ZJ/JAfyFDoWQiG1uao7uwOBsPIp3046YsXrCu/TUiF2EyfpvxYu0BkDV0hUzG4fDdLDO7fHP0BIyIWhRSD6YfZ36/3VXW0aWVAQAAl8n4/XTgNoVXKhScIQDgZWwYp6NB+e1hTVxRWJceAM9jEVmYpO75y49aiADBMu/bhVjDVab2X5/2gZCxhFLh1l7dqRO5eDZBJ+lm2hURg62jExITOnbLnKvBZJQqtluAtcaeEWQg9znT2fRGd4HpoMrj8ozIRKh22O0MJykZi0mnnxqOdOokFgS6byENbWlqraPeIGcX9JXD3hs4Iroi51kaHCWK80Lgdxs7mfG8El7o3mxcme3BtWh+fZ+MRGzARXw8hPl8ZSpocfnQu+tHhyMyJZv5/eP64C5I9YJs/SHpJs7Xq6fPj6ttSW1MOA97N59u047fvPxcrnJTueZVhkTU0ZzbqrRV5bx/nT5EZByRABxRQ39E8/QlYABABP+tfgPW3c4ebZ+CYQAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "2\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABrElEQVR4nL3S3Y/SQBAA8Jn9aCltqSBHTE9PTsldjMk9+v+/q4mY3OFHzjuiSAs97yi0lHa74wOUovHZfZnN/pKZ2cwA/P+D/3hCBproEBEQkAEQcs6xVLkGEDtiyIRpWaYhBWVZsi600ltELqS0Wt2nPQsRy1UwmS9LTSQAALlsWG73+fmpky8yaaeOSlIOBAIAkAlpP+5fvHLnP+fQ60pb6qKkbU0iQmY/8dnoarxqG8+wSJJNXqFWCp22/vr28pcwmk4Y36e5qrqlUhs2hp8uQ+32B+YivMsKggpBNIxiGeGj5uDNyV3wPc5pj8A5qCR3X5B38boxu4k2hLBHZPmywR1bds464edxrDTUCMVqurYYc176MBpOlhutkSpUa8pm3HDOe97tu1FUlAR1Q+VGpSjslt/PPnwMMg1AdVrSUDLhHJ95397fZJp2k2JbRGCy5Q+OH4ZXi73tEAl58+jk1LgeTlU9dlZdDPfI96IvtynVKKrYaHe8IriOCqK/EbllmevJ+EesCACQ/kAq01mug2lSEuw+sl8wJixTmjyJk4OSB9tHDIng0OA3lAjRT2O6MkgAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "4\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABtklEQVR4nLWSSWsbQRCFq5fZWtLMaEGWSWwRG4JP+f8/I4dA4pOxwRhZWMuMRq2e6q7KQdLYyTGQOjV8vFfF6wfwX0Z8eEoABua/oZA67g0Kk8BhV1XWExHzCSqdFrOrL9djE5rl8/3zen9oAzFoAJA6zT/f3d1OJFMxKhPtnZAEAPpoWoxHA/e4bSBOScbSY2A+QpZK4PaJVgsbj6bGtkRMLI62QGSrxQ6XS8ijnKraCgF8sgXv9tvYB5+l44u+a2rnic6QqbWp6/eGcT4xuLTeBwLoIHkM6eW0KLKwqXWkTkFoAACg4FHn0zLTBzMa72vrO1sACOjqt8hmCkNyGZxtw7uSqa1f/CqTrMzYXLXrBok7JaPdhLVs2yiffy2nw5U9xyeklFqCQ9w7WWfX5aBvYndSCqWiNDPCWduy8sIkSSTOO4WK0jzvcYM6yPJilDIF7kIQiSlmOdTZjvPbbzeuqez7tUyqP5mVaF00mc8OD4+vDXZKj9i66NM0EbqnFj++/3rd4fk/Obh6ITzifBBeHn7eP70d2alDUuu0NyiHia9Wm+qAxH+0T4AAEAIY+GP//nV+A+iU+ep185/TAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "3\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABlUlEQVR4nJWSUW8SQRSF77k7yw6w0BakdW2KTZNGH/z/f0SjNrYhjaRqaYuFwuLszhwfWLoQY6L3aWa+zD1nzh2R/ylsrc3f0e4OKgAYhELhDlRAVeO9I13MZitfcqstFKra6L45i26+gisNrKESgGkenp0HN53ZMrCGEACmtZ8Ns7mJTFKgKIJWDCKCuN3PXlj3qxC1cUTRyjIDGTUP+sOBWS49xfuyghRSJGr1uq9OWvntz/kid76s3ULjdD9Nj7rlzWh8t/Te1wkBke12Gv2TZH5x9X1RkKRUmhCoTeL07evO5Msk9wykyEZTNIqQnL7r5e/HTyWBreABFY/med9ffn4o+BzMmhkFmsMsuv9w60iuc69uakNN5+Xpnru+WIawFqzbxtYeDrPk/uOdC1zPauNWYOzB8fHAj0erIoiiGuUGNnuDzM6uJrkTCBRk3ZbGpu1yPHqYeUIk+C3NUAbRZX55PXU+BJKQ55sIbvH4Q+afvj05BmLjFuvUTZK2bTl9dP7P3weNFPS+evwuFMju8T/Ub0AawExDVNKjAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "5\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABlElEQVR4nL3SyW4TQRAG4Kre0rN5TOwktgEnyikSEuLA+z8B4oKQEBImkJCJHQ92ehbP0guHMJO5I1Gnlj79fx26AP7/4OBJqJS+IE1zqI2xxg2Qiuj09eIsdEW6ftiqrDKOdUaC0/M3VzOP4yG9/X6dWGN7RDG+eP9uaquSUi/yBGPEQY/y+PJSpnf3dfiC7lV2qE2PSI4mJ5h8XT3KmZbp9reqDXa1SIVP9a9vdyTyWL5dP5bNcxIoJdhq7yyc++1Dsiub1kJfixTAW0ZNyNXteqMq7aBDdAjOyvlSO7XtSgEIAAA4sEbnhQtiUZdZ1TxZv9OZbJc4n+oCfWasgx4RGSdtlmSUHMnRmP81YABAgArp8Xa3r/j8Ik5xiEhFEI6kzVUen8QBbzUMkjyYTEbMloXms1cqzweIRITTRUxLZMHVW391k7lBLQCV05c+5jZeju8//yzcc9JpfShaf3EsCOibDx83rR3stJXarDihEeyvv3z6oXSXxKc/8aLpZOxBniZpYTrrbogQAojWWdvLP80fg+bNIwXecvUAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "8\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABo0lEQVR4nLWSy5LaQAxFpW657TbY2IMhhMzj//8pu2QDMwYCNnY/lYUZmHWqoo0WR0fqrroA/1r4aAgCEUME5hukzyFSWZamEtxgx6ub+GSioHxRrlfrivrdrt23g+eHSbpar16evzfp8P7rp/IhPKBM53W9rFUcfLJg44ariZ8QiUgIfxo/cJ5lWNSt7v3dZI7m7E4QUS+3T6REgvfXcrTX4BQ4jvlYrRSGyI+vRD9EL2IICYRMG2td+AIDsxcQEXQ5I2+GMSLfb0ZkjoJUXjVztOfOflmLKGWiZZoUr5u5PZ46B4A8QSFTPSsLrfL19of+3V4syBgmE2Ve1M2yrIqqmun+0llJHuJkyqxYbzZNNm+qVLuT8Sr11oUbVKu31yWpnKSWhlXZeerATWtlvX1+0TJJFGGwRpUDB29lAAJAOSuXTS0gkOgv5551Cd5fEZkAgMFTtlDejPvDoQvo/Gg83G6Gy3FXJLk57nb71gtBY/fnMkYAAmDfRrP/VoyH9uPQYyIwGHsd+JYsSnSeEfbD6CIIAo7e8SNDCIwIfM/df62/z8ngRcdUFskAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "2\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABRklEQVR4nLWSS0/cQBCEq3pmMCywSLw5RuL//x6uOaBIESAeZs163F052LvrJZdc0pcZ6VOpqksN/Jfh9kMSUADQd8hkVo4vTtbtc+uaeJ5YTikvbu7v4vGhrodJO0KaWTo4Ob/9Yen30xeoPcjcLI6Xy8v8dlos9pSCkTkNXpIlbCw3npC8rt768HUdtnHzGDlcKrknEJ2L87SCQsmdR4cY2j5iUtr4REhgWV4ddW3dlmAbz4jg4uw0t52M3yCAkBrqM7Zsl5aUg+YugIQ0U5JSKt63LMVIkLNuCSmVYfURJQ8SBWhSEgBYCtarConcrw8IpOKrTiZB2vc0Qy5lCMAMJKDZKqTlxgAFbcy6WwUQi9XaD+4hBTTf0+nvL4+vP5+73mNkuwMzHhxeX+L112fVptzd9QGpaVA7//v6AACJlAv/Nn8AUVOr7xQAlfgAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "1\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<PIL.Image.Image image mode=L size=28x28>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAABwAAAAcCAAAAABXZoBIAAABuklEQVR4nLXSWZOaQBAA4J4RWMATFOLigVe5lcf8//+RLVeNRF0pWU9kBzmGzkugNpvHVPppqr/q6ZmuBvgvQT6eKRACCAiYfUIqyOW6qogk9G8BS7KPSJVaw+q3W+Xs6i5Wb36CAEJhDWvUG1gViUbbCsEogQJLD3pvPB4ZKnAuWMCYz3iORGh8GT0NdDwGAapNbXI+nMMcS4rWn4x15u39S6pPOu2+IxUoVNsj2wjX6+3uChaUTct8yHsSURtOuunrwtm+xQJog0fdLBP8jWpnMhRe57PNOeLC7fwOqiLllVS3bfX04/vyEiNiGmWAQPJKwegYfDdfHeMMgIplGaI7IAAFAJDbHcn/uTncOWZEqpvV9HKO82vlVoOeth5LCaFSszfQQ88P8/EpFSm5+AlQIkrG+Kstu67HigfxBMVKNYxFpWVPp01/5Xhh8c8g0Fo2q3BFswZ27fKydI7F4MODaxjf2tdIqD3WxcPL88x9xxyD3UqfDrsMSyLenPls4d445MjWMs2GmpIkt72zXG2O9xQKTA/86nbNKgb73dY7hTz7Y8GIJKuyhPzOwgQR8K/tIwCQp/81fgGaLNUpHaAJHgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.45800865800865803"
            ]
          },
          "metadata": {},
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oZjP6jnBgOxB"
      },
      "source": [
        "Finally, we'll rewrite the K-means model as an active learning problem and perform semi-supervised clustering of the data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Uzl-FUtPHMEx",
        "outputId": "292f3b0f-5b82-4e51-e692-e85687983aa8",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "source": [
        "from scipy import stats\n",
        "import math\n",
        "from functools import reduce\n",
        "\n",
        "# ACTIVE LEARNING VERSION\n",
        "\n",
        "#for nlabeled in range(20,len(X_train),10):\n",
        "nlabeled = 10\n",
        "if True:\n",
        "  ans = []\n",
        "  k = 50 # 10 # 2 # 5 # 20\n",
        "  if True:\n",
        "#  for k in range(10,200,50):\n",
        "    bestcenters, bestloss = kmeans(X_train, k)\n",
        "    for rep in range(9):\n",
        "      centers, loss = kmeans(X_train, k)\n",
        "      if loss < bestloss: bestcenters, bestloss = centers, loss\n",
        "    # How do we test the clustering that was discovered?\n",
        "    # Assign testing points to clusters\n",
        "    test_centerids, loss = assign_data(X_test, bestcenters)\n",
        "\n",
        "    # Let's label one example in each category\n",
        "    train_centerids, loss = assign_data(X_train, bestcenters)\n",
        "\n",
        "    clust_labs = np.repeat(labs[0],k)\n",
        "    for i in range(len(bestcenters)):\n",
        "      clust_labs[i] = y_train[train_centerids == i][0]\n",
        "\n",
        "    # semi-supervised clustering\n",
        "    print(k,sum(clust_labs[test_centerids] == y_test)/len(y_test))\n",
        "\n",
        "    # nearest neighbors\n",
        "    labids, loss = assign_data(X_test, X_train[:k])\n",
        "    print(k, sum(y_train[labids] == y_test)/len(y_test))"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "50 0.7164502164502164\n",
            "50 0.6181818181818182\n"
          ]
        }
      ]
    }
  ]
}