{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "109aea99",
   "metadata": {},
   "source": [
    "# Hands-on clustering #2:\n",
    "## Ensemble clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82577d25",
   "metadata": {},
   "source": [
    "### Imports and formatting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0545f90c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import ttest_ind\n",
    "from scipy.special import binom\n",
    "\n",
    "from sklearn.cluster import KMeans, AgglomerativeClustering, MiniBatchKMeans\n",
    "from sklearn import datasets\n",
    "from sklearn.metrics import normalized_mutual_info_score\n",
    "\n",
    "import h5py\n",
    "\n",
    "plt.style.use(\"seaborn\")\n",
    "plt.rcParams.update(\n",
    "    {\n",
    "        \"font.size\": 18,\n",
    "        \"xtick.labelsize\": 18,\n",
    "        \"ytick.labelsize\": 18,\n",
    "        \"legend.fontsize\": 20,\n",
    "        \"axes.labelsize\": 20,\n",
    "        \"axes.titlesize\": 20,\n",
    "        \"figure.figsize\": (8, 8),\n",
    "    }\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d141cd4",
   "metadata": {},
   "source": [
    "## 2. Understanding ensemble methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e408f004",
   "metadata": {},
   "source": [
    "#### Binomial distribution functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d211c30c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def P_k_np(k, n, p):\n",
    "    return binom(n, k) * p**k * (1 - p) ** (n - k)\n",
    "\n",
    "\n",
    "def P_supk_np(k, n, p):\n",
    "    P = 0\n",
    "    for i in range(k, n + 1):\n",
    "        P += P_k_np(i, n, p)\n",
    "    return P\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d13b174a",
   "metadata": {},
   "source": [
    "### 2.1 Probability of correct classification when ensembling 1001 models with an accuracy of 51%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "757f4d8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of correct prediction:  0.7366309083614334\n"
     ]
    }
   ],
   "source": [
    "p = 0.51\n",
    "n = 1001\n",
    "\n",
    "# If k models predict the correct answer, the ensemble is correct when k > n/2\n",
    "k_min = int(n / 2) + 1\n",
    "\n",
    "print(\"Probability of correct prediction: \", P_supk_np(k_min, n, p))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98fc8ec6",
   "metadata": {},
   "source": [
    "### 2.2 With 1001 models, what accuracy is needed on the weak models to outperform the expert model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1aec0784",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54320 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 1001\n",
    "k_min = int(n / 2) + 1\n",
    "\n",
    "P_ensemble = []\n",
    "x_plot = np.arange(0.50, 0.60, 0.01)\n",
    "for p in x_plot:\n",
    "    P_ensemble.append(P_supk_np(k_min, n, p))\n",
    "\n",
    "P_expert = [0.9 for _ in x_plot]\n",
    "plt.plot(x_plot, P_ensemble)\n",
    "plt.plot(x_plot, P_expert)\n",
    "plt.legend([\"Ensemble\", \"Expert\"])\n",
    "plt.xlabel(\"Accuracy of weak models\")\n",
    "plt.ylabel(\"Probability of success\")\n",
    "plt.title(\n",
    "    \"Success rate of ensemble learning with 1001 weak models for different accuracies\"\n",
    ")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2671b819",
   "metadata": {},
   "source": [
    "#### Comment\n",
    "With 1001 weak models, each with an accuracy of 0.52, ensemble performance is 0.9 (equal to expert here).\n",
    "\n",
    "Hence, with an accuracy > 0.52, the ensemble outperforms the expert."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f071e0d1",
   "metadata": {},
   "source": [
    "### 3.3 To work as expected, three conditions must be verified by the set o weak learners"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "131310e7",
   "metadata": {},
   "source": [
    "#### 1. The number of weak learners must be \"large\"\n",
    "\n",
    "The following example presents the ensemble accuracy for different number of weak learners when the individual accuracy is 0.53"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3aebc1bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 58004 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = 0.53\n",
    "\n",
    "P_ensemble = []\n",
    "x_plot = np.arange(51, 1011, 10)\n",
    "for n in x_plot:\n",
    "    k_min = int(n / 2) + 1\n",
    "    P_ensemble.append(P_supk_np(k_min, n, p))\n",
    "\n",
    "P_expert = [0.9 for _ in x_plot]\n",
    "plt.plot(x_plot, P_ensemble)\n",
    "plt.plot(x_plot, P_expert)\n",
    "plt.legend([\"Ensemble\", \"Expert\"])\n",
    "plt.xlabel(\"Number of weak models\")\n",
    "plt.ylabel(\"Probability of success\")\n",
    "plt.title(\n",
    "    \"Success rate of ensemble learning with 0.53 accuracies for different number of learners\"\n",
    ")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9ad01f6",
   "metadata": {},
   "source": [
    "#### 2. The individual learners must be \"good\"\n",
    "\n",
    "Individual models in ensembles must be at least better than random classification.\n",
    "\n",
    "The following example illustrates the behavior of ensemble when individual models have an accuracy of 0.47"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "fc9b0f2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 51456 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = 0.47\n",
    "\n",
    "P_ensemble = []\n",
    "x_plot = np.arange(51, 1011, 10)\n",
    "for n in x_plot:\n",
    "    P_ensemble.append(P_supk_np(int(n / 2) + 1, n, p))\n",
    "\n",
    "P_expert = [0.9 for _ in x_plot]\n",
    "plt.plot(x_plot, P_ensemble)\n",
    "plt.plot(x_plot, P_expert)\n",
    "plt.legend([\"Ensemble\", \"Expert\"])\n",
    "plt.xlabel(\"Number of weak models\")\n",
    "plt.ylabel(\"Probability of success\")\n",
    "plt.title(\n",
    "    \"Success rate of ensemble learning with 0.53 accuracies for different number of learners\"\n",
    ")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e34cd80",
   "metadata": {},
   "source": [
    "Obviously, there is a strong relation between how good the individual learners are and the number of weak models that is needed for good performance. For example, if weak models have an accuracy of 0.7, 13 models are sufficient to generate strong ensemble predictions (> 0.9)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ac994535",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of correct prediction:  0.9376247882008\n"
     ]
    }
   ],
   "source": [
    "p = 0.7\n",
    "n = 13\n",
    "k_min = int(n / 2) + 1\n",
    "\n",
    "print(\"Probability of correct prediction: \", P_supk_np(k_min, n, p))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f462f1e4",
   "metadata": {},
   "source": [
    "#### 3. Weak learners must be diverse: \n",
    "\n",
    "If all weak model predict the same thing, i.e., make the same mistakes, then the ensemble is not better than the individual weak models. \n",
    "\n",
    "_Note: All the computations above assume that all weak models are independent, which is rarely true in practice_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a84b2325",
   "metadata": {},
   "source": [
    "## 3. Ensemble methods for clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d4cc40c",
   "metadata": {},
   "source": [
    "### 3.1 Generating clustering partitions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b69e6eff",
   "metadata": {},
   "source": [
    "#### List of possible ways to generate clustering partitions\n",
    "\n",
    "1. Use different clustering algorithms\n",
    "    1. Different approaches\n",
    "    2. Different parameters\n",
    "2. Use different features\n",
    "    1. Subsets of the features\n",
    "    2. Different modeling of the problem\n",
    "    3. Different data transformation\n",
    "3. Use different initializations\n",
    "4. Use different subsets of data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afe8bfa0",
   "metadata": {},
   "source": [
    "#### Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fe7e42d9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 83468 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_samples = 500\n",
    "noisy_circles, labels_noisy_circles = datasets.make_circles(\n",
    "    n_samples=n_samples, factor=0.6, noise=0.05\n",
    ")\n",
    "\n",
    "plt.scatter(\n",
    "    noisy_circles[:, 0], noisy_circles[:, 1], c=labels_noisy_circles, cmap=\"tab20\"\n",
    ")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c798ca8",
   "metadata": {},
   "source": [
    "#### 3.1.1 Run regular Kmeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e0b9d09c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 81276 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NMI:  4.6166733759236576e-05\n"
     ]
    }
   ],
   "source": [
    "kmeans = KMeans(n_clusters=2)\n",
    "clusters = kmeans.fit_predict(noisy_circles)\n",
    "\n",
    "plt.scatter(noisy_circles[:, 0], noisy_circles[:, 1], c=clusters, cmap=\"tab20\")\n",
    "plt.show()\n",
    "\n",
    "print(\"NMI: \", normalized_mutual_info_score(labels_noisy_circles, clusters))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e29b727",
   "metadata": {},
   "source": [
    "#### Comments:\n",
    "\n",
    "Kmeans clearly did not work for these data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ff4ecf2",
   "metadata": {},
   "source": [
    "#### 3.1.2 Generating 100 partitions using Kmeans with 100 clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8438ed79",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=100)\n",
    "\n",
    "clusters = []\n",
    "for _ in range(100):\n",
    "    clusters.append(kmeans.fit_predict(noisy_circles))\n",
    "clusters = np.array(clusters)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1a97583",
   "metadata": {},
   "source": [
    "Let's now display 3 of these cluster assignments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2696b7e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 233372 characters]",
      "text/plain": [
       "<Figure size 1728x576 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(24, 8))\n",
    "for i in range(3):\n",
    "    plt.subplot(1, 3, i + 1)\n",
    "    plt.scatter(noisy_circles[:, 0], noisy_circles[:, 1], c=clusters[i], cmap=\"tab20\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93f33467",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "- Clusters are definitely \"good\" (better than random)\n",
    "- There is \"a high number\" of partitions (> 100)\n",
    "- Regarding diversity, all three clustering partitions shown above are slightly different but most points appear to share common clustering within the three partitions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bd5d602",
   "metadata": {},
   "source": [
    "### 3.2 Consensus function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8dfcc2",
   "metadata": {},
   "source": [
    "- The grouping strategy presented in Section 2 (simple majority voting) is used for supervised classification problems. For clustering, it cannot be used straightforwardly because of the label correspondance problem.\n",
    "- The most straightforward approach to consensus would be relabelling and voting. Other approaches exist, here we will see one based on counting co-occurrence of points in clusters between different partitions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d9fde89",
   "metadata": {},
   "source": [
    "#### 3.2.1 Co-Association Matrix function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ff6a6ed0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_CoAssociationMatrix(clusters):\n",
    "    n_samples = clusters.shape[1]\n",
    "    CAM = np.zeros((n_samples, n_samples))\n",
    "    for i in range(clusters.shape[0]):\n",
    "        for j in range(n_samples):\n",
    "            for k in range(j + 1, n_samples):\n",
    "                if clusters[i][j] == clusters[i][k]:\n",
    "                    CAM[j, k] += 1\n",
    "                    CAM[k, j] += 1\n",
    "    return CAM\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "219e25a7",
   "metadata": {},
   "source": [
    "#### 3.2.2 Cluster the Co-Association Matrix\n",
    "\n",
    "- Agglomerative clustering can be used to cluster a dataset based on precomputed distances\n",
    "- The Co-Association Matrix represents the similarity between two points in dataset. As Agglomerative Clustering requires distances (dissimilarities) we need to pass the opposite of the Co-Association Matrix as argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "abeb4f13",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 83468 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "CoAssociationMatrix = compute_CoAssociationMatrix(clusters)\n",
    "\n",
    "agg = AgglomerativeClustering(n_clusters=2, affinity=\"precomputed\", linkage=\"average\")\n",
    "clusters = agg.fit_predict(-CoAssociationMatrix)\n",
    "\n",
    "plt.scatter(noisy_circles[:, 0], noisy_circles[:, 1], c=clusters, cmap=\"tab20\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0dd2296",
   "metadata": {},
   "source": [
    "#### 3.2.3 Initialization issue\n",
    "\n",
    "- The algorithm did not work as expected\n",
    "- Previously, we noted that the diversity in the partitions was limited, this is due to the default parameters of Kmeans. Indeed, by default, Scikit-learn Kmeans runs the \"Kmeans++\" algorithm to properly initialize parameters and conducts multiple initializations (10). This has the effect of limiting the diversity in the partitions because each run converges to a similar local minimum.\n",
    "- To correct this issue we need to set Kmeans initialization to random and the number of restart to 1. this will \"reduce the performance\" of the individual partitions, but increase the diversity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3c700bb1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 83468 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kmeans = KMeans(n_clusters=100, init=\"random\", n_init=1)\n",
    "\n",
    "clusters = []\n",
    "for _ in range(100):\n",
    "    clusters.append(kmeans.fit_predict(noisy_circles))\n",
    "clusters = np.array(clusters)\n",
    "\n",
    "CoAssociationMatrix = compute_CoAssociationMatrix(clusters)\n",
    "\n",
    "agg = AgglomerativeClustering(n_clusters=2, affinity=\"precomputed\", linkage=\"average\")\n",
    "clusters = agg.fit_predict(-CoAssociationMatrix)\n",
    "\n",
    "plt.scatter(noisy_circles[:, 0], noisy_circles[:, 1], c=clusters, cmap=\"tab20\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35af2f72",
   "metadata": {},
   "source": [
    "#### 3.2.4 Find the good number of clusters in each partitions for 100 partitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "4207b096",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 NMI:  0.00011566251563091068\n",
      "10 NMI:  0.00023373198039611758\n",
      "12 NMI:  0.0\n",
      "14 NMI:  0.00214617865119764\n",
      "16 NMI:  0.0024818544496891283\n",
      "18 NMI:  0.006031990097249684\n",
      "20 NMI:  0.9999999999999974\n",
      "22 NMI:  0.9999999999999974\n",
      "24 NMI:  0.9999999999999974\n"
     ]
    }
   ],
   "source": [
    "nmis = []\n",
    "for k in range(8, 26, 2):\n",
    "    kmeans = KMeans(n_clusters=k, init=\"random\", n_init=1)\n",
    "\n",
    "    clusters = []\n",
    "    for _ in range(100):\n",
    "        clusters.append(kmeans.fit_predict(noisy_circles))\n",
    "    clusters = np.array(clusters)\n",
    "\n",
    "    CoAssociationMatrix = compute_CoAssociationMatrix(clusters)\n",
    "\n",
    "    agg = AgglomerativeClustering(\n",
    "        n_clusters=2, affinity=\"precomputed\", linkage=\"average\"\n",
    "    )\n",
    "    clusters = agg.fit_predict(-CoAssociationMatrix)\n",
    "\n",
    "    nmis.append(normalized_mutual_info_score(labels_noisy_circles, clusters))\n",
    "\n",
    "    print(k, \"NMI: \", nmis[-1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5d493dd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fca08d6e9d0>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 17960 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(8, 26, 2), nmis)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ddc626d",
   "metadata": {},
   "source": [
    "#### Comments:\n",
    "Ensemble of 100 partitions containing 20 clusters appear to work well for the concentric circle data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ab5df10",
   "metadata": {},
   "source": [
    "### 3.3 Repeat the process for the elongated blob data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "41ef8c72",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fca0955e6a0>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 54948 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "random_state = 170\n",
    "X, labels_blobs = datasets.make_blobs(n_samples=n_samples, random_state=random_state)\n",
    "transformation = [[0.6, -0.6], [-0.4, 0.8]]\n",
    "blobs = np.dot(X, transformation)\n",
    "\n",
    "plt.scatter(blobs[:, 0], blobs[:, 1], c=labels_blobs, cmap=\"tab20\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6a9c5ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=20, init=\"random\", n_init=1)\n",
    "\n",
    "clusters = []\n",
    "for _ in range(100):\n",
    "    clusters.append(kmeans.fit_predict(blobs))\n",
    "clusters = np.array(clusters)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e0dda8f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54548 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "CoAssociationMatrix = compute_CoAssociationMatrix(clusters)\n",
    "\n",
    "agg = AgglomerativeClustering(n_clusters=3, affinity=\"precomputed\", linkage=\"average\")\n",
    "clusters = agg.fit_predict(-CoAssociationMatrix)\n",
    "\n",
    "plt.scatter(blobs[:, 0], blobs[:, 1], c=clusters, cmap=\"tab20\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccd385ae",
   "metadata": {},
   "source": [
    "## 4. Practical implementation of Ensemble clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "102ff7bb",
   "metadata": {},
   "source": [
    "- The objective is to group together images representing the same person. The dataset contains 2800 images (200 persons and 14 images per person). As there is 200 persons, the number of expected clusters is 200."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16c718c0",
   "metadata": {},
   "source": [
    "#### 4.1 Running Kmeans on a single representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "97ed54f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"representations/densenet121.h5\"\n",
    "\n",
    "f = h5py.File(filename, \"r\")\n",
    "data = np.array(f[\"data\"])\n",
    "labels = np.array(f[\"labels\"])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "79acab21",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NMI:  0.8851436238144947\n"
     ]
    }
   ],
   "source": [
    "km = MiniBatchKMeans(200)\n",
    "clusters = km.fit_predict(data)\n",
    "\n",
    "print(\"NMI: \", normalized_mutual_info_score(labels, clusters))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84971b19",
   "metadata": {},
   "source": [
    "#### 4.2 Running Kmeans on all representations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f09a04e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "directory = \"representations/\"\n",
    "files = [\n",
    "    \"densenet121.h5\",\n",
    "    \"densenet169.h5\",\n",
    "    \"densenet201.h5\",\n",
    "    \"inception.h5\",\n",
    "    \"inception_resnet.h5\",\n",
    "    \"nasnet.h5\",\n",
    "    \"resnet.h5\",\n",
    "    \"vgg16.h5\",\n",
    "    \"vgg19.h5\",\n",
    "    \"xception.h5\",\n",
    "]\n",
    "\n",
    "km = MiniBatchKMeans(200)\n",
    "clusters = []\n",
    "for f in files:\n",
    "    ff = h5py.File(directory + f, \"r\")\n",
    "    data = np.array(ff[\"data\"])\n",
    "    labels = np.array(ff[\"labels\"])\n",
    "    clusters.append(km.fit_predict(data))\n",
    "clusters = np.array(clusters)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6cba086a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "densenet121.h5\n",
      "        NMI:  0.8817267641060259\n",
      "densenet169.h5\n",
      "        NMI:  0.8788637674639859\n",
      "densenet201.h5\n",
      "        NMI:  0.8794119764315856\n",
      "inception.h5\n",
      "        NMI:  0.8952033782797599\n",
      "inception_resnet.h5\n",
      "        NMI:  0.8385824811497752\n",
      "nasnet.h5\n",
      "        NMI:  0.8866301868250733\n",
      "resnet.h5\n",
      "        NMI:  0.8700830295425884\n",
      "vgg16.h5\n",
      "        NMI:  0.8444703426760213\n",
      "vgg19.h5\n",
      "        NMI:  0.8845026162210827\n",
      "xception.h5\n",
      "        NMI:  0.8718054429574429\n"
     ]
    }
   ],
   "source": [
    "for i, f in enumerate(files):\n",
    "    print(f)\n",
    "    print(\"        NMI: \", normalized_mutual_info_score(labels, clusters[i]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70e0899a",
   "metadata": {},
   "source": [
    "#### Ensemble the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "950c9cf9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ensemble\n",
      "        NMI:  0.9372664229280685\n"
     ]
    }
   ],
   "source": [
    "CoAssociationMatrix = compute_CoAssociationMatrix(clusters)\n",
    "\n",
    "agg = AgglomerativeClustering(n_clusters=200, affinity=\"precomputed\", linkage=\"average\")\n",
    "clusters_ensemble = agg.fit_predict(-CoAssociationMatrix)\n",
    "\n",
    "print(\"ensemble\")\n",
    "print(\"        NMI: \", normalized_mutual_info_score(labels, clusters_ensemble))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "afe6880b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 34920 characters]",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_bars = [i + 1 for i in range(len(clusters))]\n",
    "y_bars = [normalized_mutual_info_score(labels, c) for c in clusters]\n",
    "\n",
    "plt.bar(x=x_bars, height=y_bars)\n",
    "plt.bar(x=[11], height=normalized_mutual_info_score(labels, clusters_ensemble))\n",
    "plt.ylim([0.85, 0.95])\n",
    "plt.xticks([i + 1 for i in range(11)], labels=files + [\"Ensemble\"], rotation=90)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1b7e9e8",
   "metadata": {},
   "source": [
    "## 5. Going further"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a257f41",
   "metadata": {},
   "source": [
    "Here, we propose to test combining"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "4a5ffd91",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(150, 4)\n",
      "(150,) {0, 1, 2}\n"
     ]
    }
   ],
   "source": [
    "data = datasets.load_iris()\n",
    "X = data.data\n",
    "y = data.target\n",
    "\n",
    "print(X.shape)\n",
    "print(y.shape, set(y))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "6c2f0b42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NMI kmeans:  0.7581756800057782\n",
      "NMI agglomerative:  0.7700836616487866\n",
      "NMI spectral:  0.7979885217013317\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cluster import DBSCAN, SpectralClustering, AffinityPropagation, MeanShift\n",
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "kmeans = KMeans(n_clusters=3)\n",
    "clusters_km = kmeans.fit_predict(X)\n",
    "print(\"NMI kmeans: \", normalized_mutual_info_score(y, clusters_km))\n",
    "\n",
    "agg = AgglomerativeClustering(n_clusters=3)\n",
    "clusters_agg = agg.fit_predict(X)\n",
    "print(\"NMI agglomerative: \", normalized_mutual_info_score(y, clusters_agg))\n",
    "\n",
    "spectral = SpectralClustering(n_clusters=3)\n",
    "clusters_spectral = spectral.fit_predict(X)\n",
    "print(\"NMI spectral: \", normalized_mutual_info_score(y, clusters_spectral))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5bbf35f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NMI:  0.7979885217013317\n"
     ]
    }
   ],
   "source": [
    "clusters = [\n",
    "    clusters_km,\n",
    "    clusters_agg,\n",
    "    clusters_spectral,\n",
    "]  # , clusters_affinity, clusters_shift, clusters_db, clusters_mixt]\n",
    "clusters = np.array(clusters)\n",
    "\n",
    "CoAssociationMatrix = compute_CoAssociationMatrix(clusters)\n",
    "\n",
    "agg = AgglomerativeClustering(n_clusters=3, affinity=\"precomputed\", linkage=\"average\")\n",
    "clusters_ensemble = agg.fit_predict(-CoAssociationMatrix)\n",
    "\n",
    "print(\"NMI: \", normalized_mutual_info_score(y, clusters_ensemble))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10d768ca",
   "metadata": {},
   "source": [
    "#### Comments:\n",
    "Although ensemble did not allow to reach a better classification, it still allowed to reach a partition corresponding to the best algorithm (here spectral), which was not known a priori."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e24cf64",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
