transfer_learning/kmeans_CNN_clustering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# kmeans_CNN_clustering\n",
    "- In [CNN_transfer_learning](./CNN_transfer_learning), we used a pre-trained convolutional neural network to generate a new feature set to better represent the galaxy images. These features were then compressed by a principal component analysis.\n",
    "- Here, we take the output from the PCA to cluster together galaxy images using the K-means clustering algorithm.\n",
    "- We select 10 clusters as a baseline, before exploring more complex clustering algorithms. \n",
    "- Using output from CNN with 128 vs 37632 output features (compressed to 20 vs 100). Both return sensible clustering, potentially better performance with more parameters (however needs to be quantified)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np \n",
    "import seaborn as sns\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "import cluster_evaluate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading set of features \n",
    "- Also weights in case we want to manipulate feature space to rank pca components that best describe the variance in the data\n",
    "- Several feature representations were generated, so selecting combination of output features from CNN and the number of features these were compressed to by the pca."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'./cluster_plots/kmeans/CNN-37632_PCA-100/'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca_feat = 100\n",
    "cnn_feat = 37632\n",
    "\n",
    "features = np.load('./features/CNN-' + str(cnn_feat) + '_PCA-' + str(pca_feat) + '_features.npy')\n",
    "weights = np.load('./features/CNN-' + str(cnn_feat) + '_PCA-' + str(pca_feat) + '_weights.npy')\n",
    "\n",
    "plot_path = './cluster_plots/kmeans/CNN-' + str(cnn_feat) + '_PCA-' + str(pca_feat) +'/'\n",
    "plot_path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating K-means cluster object and fitting to the set of features output from the pca\n",
    "- Generating count plot of cluster sizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KMeans(n_clusters=10)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kmeans = KMeans(n_clusters=10)\n",
    "kmeans.fit(features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEKCAYAAADn+anLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXZklEQVR4nO3dfbRddX3n8fdHAuVBEZXAYAIGK0VcCIJZSGGJSmoVRaAWLFoRKQ7WgoOPFeh0xKkua4tVGJ+KoMLIiIi4YFpHRRGfRRNAASMSUSCAEgQfADEg3/lj7yuXcHdyQs7Z53Dzfq111j374Zzfl6yQz92/396/X6oKSZJm8ohxFyBJmlyGhCSpkyEhSepkSEiSOhkSkqROc8ZdwDBtueWWtWDBgnGXIUkPK0uWLLm1qubOdGxWhcSCBQtYvHjxuMuQpIeVJNd1HbO7SZLUyZCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktRpVj1xLa3J+974f3tr65h3v6i3tqRR8UpCktTJkJAkdbK7SRqDd7z84N7a+oePn9tbW5p9vJKQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ1m7S2wT3/zmb21teRfX9FbW5LUJ68kJEmdeg2JJK9PclWSK5N8IsnGSbZPckmSa5J8MslG7bl/1G4va48v6LNWSVKPIZFkHvDfgIVVtTOwAXAo8C7gPVW1A3A7cGT7kSOB26vqScB72vMkST3qu7tpDrBJkjnApsDNwL7A1LwBZwAHte8PbLdpjy9Kkh5rlaT1Xm8hUVU3AicB19OEw6+AJcAvq+re9rTlwLz2/Tzghvaz97bnP27V701yVJLFSRavWLFitP8RkrSe6bO76TE0VwfbA48HNgP2m+HUmvrIao7dv6Pq1KpaWFUL586dO6xyJUn02930Z8BPqmpFVd0DnAfsBWzRdj8BzAduat8vB7YFaI8/Gritx3olab3XZ0hcD+yZZNN2bGER8APgy8DUvMmHA+e37y9ot2mPX1RVD7qSkCSNTp9jEpfQDEBfClzRtn0q8BbgDUmW0Yw5nN5+5HTgce3+NwDH9VWrJKnR6xPXVfVW4K2r7L4W2GOGc+8GDumjLknSzHziWpLUyZCQJHWatRP8SXp4OPHEE2dlW7OFVxKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROA4VEku2SZIb9SbLd8MuSJE2CQa8kfgLMnWH/Y9tjkqRZaNCQCFAz7H8kcPfwypEkTZI5qzuY5JT2bQHvTHLXtMMbAHsAl4+oNknSmK02JICntj8D7ASsnHZsJXApcNII6pIkTYDVhkRVPQcgyUeBY6vq171UJUmaCGu6kgCgqo4YdSGSpMkzUEgk2Rg4FlgEbMUqA95VtcuA37MFcBqwM804x98AVwOfBBYAPwVeUlW3t7fcngy8ALgLeGVVXTpIO5Kk4RgoJIAPAH8BfAr4JjPf6TSIk4HPVdXBSTYCNgVOAL5UVf+c5DjgOOAtwH7ADu3rGcAH25+SpJ4MGhIHAYdU1RcfakNJNgf2AV4JUFUrgZVJDgSe3Z52BnAxTUgcCJxZVQV8O8kWSbapqpsfag2SpLUz6HMSdwE3rGNbTwRWAB9NclmS05JsBmw99Q9/+3Or9vx5q7S5vN33AEmOSrI4yeIVK1asY4mSpOkGvZL4F+ANSV5TVfetQ1u7A6+tqkuSnEzTtdTlQdOAMEM3V1WdCpwKsHDhwofaDTYy1//Pp675pCHZ7n9c0VtbD8VX9nlWb20966tf6a0taTYbNCSeCzwTeH6SHwD3TD9YVQcM8B3LgeVVdUm7fS5NSPx8qhspyTbALdPO33ba5+cDNw1YryRpCAbtbroV+AxwEfAz4BervNaoqn4G3JBkx3bXIuAHwAXA4e2+w4Hz2/cXAK9oJxHcE/iV4xGS1K++n5N4LXBWe2fTtcARNEF1TpIjgeuBQ9pzP0tz++symjERn9WQpJ4N2t00FFV1ObBwhkOLZji3gKNHXpQkqdOgD9NdwWqejRj0YTpJ0sPLoFcS566yvSHwNGBv4P1DrUiSNDEGHZN420z7k7wZeMJQK5IkTYx1XeP6POCvh1GIJGnyrGtI7ENz55EkaRYadOD6glV3AdsAuwEzdkVJkh7+Bh24XvWBufuAq4ATquoLwy1JkjQpXHRIWk8tfcdFvbW10z/s21tbGq61epguyROBp9A8M7G0qq4dSVWSpIkw6JjE5sDpwF/SdDW1u/Np4Miq+s2I6pMkjdGgdzedDOwCPAfYpH0tave9dzSlSZLGbdCQOAB4VVV9paruaV8XA0fRrFonSZqFBg2JTZh5SvDbgI2HV44kaZIMGhLfAP4pyaZTO9qlR98GfHMUhUmSxm/Qu5teD3wOuDHJ92nubtoVuBN43ohq05Ds/b/27q2tb7z2G721JWn0Bn1O4sokOwAvB55M88T1x4Gzquq3I6xPkjRGAz8n0YbBh0dYiyRpwgw0JpHkHUn+dob9f5vkn4ZfliRpEgw6cH0YcNkM+5cArxheOZKkSTJoSGwFrJhh/y+ArYdXjiRpkgwaEtcDz5xh/z7A8uGVI0maJIMOXP878J4kGwFTU0cuAt4JvGsUhUmSxm/QW2DfnWRL4BRgo3b3SuDkqvqXURUnSRqvtbkF9vgkb6eZKjzAD6rqjpFVJkkau7VaT6Kq7gS+O6JaJEkTZtCBa0nSesiQkCR1WqvuJkmajc751B69tfWSQ77TW1vDMOi0HIaJJK2HBu1uujnJSUl2Gmk1kqSJMmhInADsBVyZ5FtJjkzyyBHWJUmaAAOFRFV9uKr2AnYGvg68nebq4iNJ+lvRRpLUq7W6u6mqllbVm4H5NFcXLwO+muSH7bTh3i0lSbPIWg1It3M3vRj4G2BfmquK04HHA/8IPBs4dLglSpLGZaCQSLI7TTC8FLgHOBM4uqqumXbOl4CvjaJISdJ4DHol8R3gQuAo4PyquneGc5YCZw+rMEnS+A0aEn9cVdet7oR2Xqcj1r0kSdKkGHSg+eldB5K8ZUi1SJImzKAhcVaS05JsOrUjyfwkFwOvH0llkqSxGzQkngHsCVyeZGGSvwKuAH4L7Lo2DSbZIMllSf6j3d4+ySVJrknyyfYOKpL8Ubu9rD2+YG3akSStu0FXpvt+koXAB4BvAQW8qapOeQhtHkszyL15u/0u4D1VdXaSDwFHAh9sf95eVU9Kcmh73l89hPYk6WFh13M/31tb3zv4eQOdtzYPv+0KPAtYRrN06R5JHrU2RSWZD7wQOK3dDs3zFue2p5wBHNS+P7Ddpj2+qD1fktSTQWeB/Ufgq8D5NGHxdGBH4Iokz1yL9t4L/D1wX7v9OOCX026pXQ7Ma9/PA24AaI//qj1/1dqOSrI4yeIVK1asRSmSpDUZ9EriNcCLquoNVbWyqq4G/pTmuYgvDvIFSfYHbqmqJdN3z3BqDXDs/h1Vp1bVwqpaOHfu3EFKkSQNaNDnJHapqlun72h/uz8uyWcH/I69gQOSvADYmGZM4r3AFknmtN83H7ipPX85sC2wvF3P4tHAbQO2JUkagkFngb11Nce+OuB3HF9V86tqAc38ThdV1V8DXwYObk87nKZLC+CCdpv2+EVV9aArCUnS6Aw8wV+S59DM3bQdsNH0Y1W17zrU8Bbg7CRvBy6jmTCQ9uf/TrKM5grCiQMlqWeDTvD3SuBDwGdoZno9H/gTYHvg42vbaFVdDFzcvr8WeNACs1V1N3DI2n63JGl4Bh24fhNwTFVNzQJ7fFXtRhMQd4yqOEnSeA0aEk/k/ruYfgdMLV36PuCVQ65JkjQhBg2JXwBTD87dSLOMKTTPLWwy7KIkSZNh0IHrrwF/TjNf0znAKUmeCyyiWWdCkjQLDRoSx9A82wDwTuBemucezgHePoK6JEkTYNAJ/m6b9v4+msn2JEmzXGdIJHnsoF8yPUQkSbPH6q4kbmWGuZJWkfacDYZWkSRpYqwuJJ7TWxWSpInUGRJV9ZU+C5EkTZ6B526akuS/8OC5m64fWkWSpIkx6NxNjwZOAV7CKgHRckxCkmahQZ+4PolmRbqDgLuBlwFvplnzwXWnJWmWGrS7aT/gpVX1tSS/B5ZU1SeT3Ay8mvvXqJYkzSKDXklsAVzXvp++1vS3gL2GXZQkaTIMGhI/ppkJFmApcGiSAC/GJUUladYaNCQ+BuzSvv9nmi6mlcC/4hQdkjRrDTp303umvb8oyZOBhcA1VXXFqIqTJI3XWj8nAX94LsJnIyRpllttSCTZDHhRVZ3dbn+A+6cMB/g98LqqunN0JUqSxmVNYxJH0DxAN+Uw4AnA3Pb1PODo0ZQmSRq3NYXEocCZq+z7r1X1oqp6EfAW4C9HUpkkaezWFBI7AD+atv1Lmi6mKYuBnYZdlCRpMqxp4Hpz4L6pjaradobPbzjsoiRJk2FNVxI3AE9dzfFd23MkSbPQmkLiP4ETk2y86oH2zqe3tudIkmahNXU3vZPm7qark7yP+8cnngwcQxMy7xxdeZKkcVptSFTVLUn2Aj5EMx1Hpg4BXwD+rqpuGW2JkqRxWeMT11V1HbBfksfQ3O0EsKyqnNhPkma5gaflqKrbge+MsBZJ0oQZdBZYSdJ6yJCQJHUyJCRJnQwJSVInQ0KS1MmQkCR1MiQkSZ0MCUlSJ0NCktSpt5BIsm2SLydZmuSqJMe2+x+b5MIk17Q/H9PuT5JTkixL8v0ku/dVqySp0eeVxL3AG6tqJ2BP4OgkTwGOA75UVTsAX2q3AfajmStqB+Ao4IM91ipJoseQqKqbq+rS9v1vgKXAPOBA4Iz2tDOAg9r3BwJnVuPbwBZJtumrXknSmMYkkiwAdgMuAbauqpuhCRJgq/a0eTxw1bvl7b5Vv+uoJIuTLF6xYsUoy5ak9U7vIZHkkcCngddV1a9Xd+oM++pBO6pOraqFVbVw7ty5wypTkkTPIZFkQ5qAOKuqzmt3/3yqG6n9ObWI0XJg22kfnw/c1FetkqR+724KcDqwtKr+bdqhC4DD2/eHA+dP2/+K9i6nPYFfTXVLSZL6MfCiQ0OwN3AYcEWSy9t9J9Asi3pOkiOB64FD2mOfBV4ALAPuAo7osVZJEj2GRFV9nZnHGQAWzXB+AUePtChJ0mr5xLUkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjoZEpKkToaEJKmTISFJ6mRISJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnqZEhIkjpNdEgkeX6Sq5MsS3LcuOuRpPXNxIZEkg2A9wP7AU8BXprkKeOtSpLWLxMbEsAewLKquraqVgJnAweOuSZJWq+kqsZdw4ySHAw8v6pe1W4fBjyjqo5Z5byjgKPazR2Bq9ex6S2BW9fxO9bVJNQAk1HHJNQAk1HHJNQAk1HHJNQAk1HHMGp4QlXNnenAnHX84lHKDPselGhVdSpw6tAaTRZX1cJhfd/DtYZJqWMSapiUOiahhkmpYxJqmJQ6Rl3DJHc3LQe2nbY9H7hpTLVI0nppkkPiu8AOSbZPshFwKHDBmGuSpPXKxHY3VdW9SY4BPg9sAHykqq7qoemhdV2tg0moASajjkmoASajjkmoASajjkmoASajjpHWMLED15Kk8Zvk7iZJ0pgZEpKkToZEaxKmAEnykSS3JLlyHO23NWyb5MtJlia5KsmxY6pj4yTfSfK9to63jaOOtpYNklyW5D/GWMNPk1yR5PIki8dUwxZJzk3yw/bvx5+OoYYd2z+Dqdevk7xuDHW8vv17eWWSTyTZuO8a2jqObWu4alR/Do5J8IcpQH4EPJfm1tvvAi+tqh/0XMc+wB3AmVW1c59tT6thG2Cbqro0yaOAJcBBY/izCLBZVd2RZEPg68CxVfXtPutoa3kDsBDYvKr277v9toafAguramwPbiU5A/haVZ3W3nG4aVX9coz1bADcSPOQ7XU9tjuP5u/jU6rqt0nOAT5bVR/rq4a2jp1pZqLYA1gJfA54TVVdM8x2vJJoTMQUIFX1VeC2vttdpYabq+rS9v1vgKXAvDHUUVV1R7u5Yfvq/TeaJPOBFwKn9d32JEmyObAPcDpAVa0cZ0C0FgE/7jMgppkDbJJkDrAp43mGayfg21V1V1XdC3wF+IthN2JINOYBN0zbXs4Y/mGcNEkWALsBl4yp/Q2SXA7cAlxYVeOo473A3wP3jaHt6Qr4QpIl7VQ0fXsisAL4aNv1dlqSzcZQx3SHAp/ou9GquhE4CbgeuBn4VVV9oe86gCuBfZI8LsmmwAt44APIQ2FINAaaAmR9kuSRwKeB11XVr8dRQ1X9vqqeRvO0/R7t5XVvkuwP3FJVS/pst8PeVbU7zazIR7ddk32aA+wOfLCqdgPuBMY2fX/b3XUA8KkxtP0Ymp6G7YHHA5sleXnfdVTVUuBdwIU0XU3fA+4ddjuGRMMpQKZpxwA+DZxVVeeNu562W+Ni4Pk9N703cEA7HnA2sG+Sj/dcAwBVdVP78xbgMzRdpH1aDiyfdjV3Lk1ojMt+wKVV9fMxtP1nwE+qakVV3QOcB+w1hjqoqtOraveq2oemq3qo4xFgSExxCpBWO2B8OrC0qv5tjHXMTbJF+34Tmv8xf9hnDVV1fFXNr6oFNH8nLqqq3n9jTLJZexMBbRfPn9N0NfSmqn4G3JBkx3bXIqDXmxlW8VLG0NXUuh7YM8mm7f8vi2jG7nqXZKv253bAixnBn8nETsvRpzFOAfIAST4BPBvYMsly4K1VdXrPZewNHAZc0Y4HAJxQVZ/tuY5tgDPaO1geAZxTVWO7BXXMtgY+0/x7xBzg/1TV58ZQx2uBs9pfpK4FjhhDDbT9788FXj2O9qvqkiTnApfSdO9cxvim5/h0kscB9wBHV9Xtw27AW2AlSZ3sbpIkdTIkJEmdDAlJUidDQpLUyZCQJHUyJKS10M7G+qZx1yH1xZCQpkmydZKTk/w4ye+S3Jjk/yV5wYja+9g4pyCX1sSH6aRWO6HhN4DfAMfTzIXzCJonaj8EbDeu2taknY309+WDTxoyrySk+32AZrLHhVV1TlVdXVVLq+p9wK4zfSBJJTl4lX0P6JJK8uokP0pyd5IVST6fZE6SE4HDgRe231NJnt1+Zl6Ss5Pc3r7+M8kO077zxHaxmVcm+THwO2Dcs7JqFvJKQgKSPJZmAsH/Pm0diz94qNMdJFkIvJ8mDL4ObAHs2x4+iWZNgMfSTIUCcFs77cSXgW8Cz6JZUOZNwBeT7FRVd7Xnbg+8DDikPefuh1KjtDqGhNR4Es1VxLAnatuOZlrtC9pFnK6j6cYCuCPJb4HftRPoAdBOOx3giKnuoySvpllXY3/gnPbUjYDDxjQTqtYThoTUmGlNkWG4kCYYfpLk88AXgPPawOjydJqrhN+0k/pN2RT442nbyw0IjZpjElLjGpqFpnZay88VDw6YDf9wsAmD3YGX0EwxfTzwwySPX813PgK4HHjaKq8/Af592nl3rmWt0lozJCSgqm6jmSr+mHZVvgeYWttiBitopjWfOm/r6dvtd99bVRdV1fHALjQDzPu3h1fSTE8/3aU03V+3VtWyVV5jXQNd6x9DQrrf39FcFSxOckiSHZM8OclrgO93fOYimuVEFybZDfgY0waQk+yf5NgkuyV5As1A86O4f+zjp8DObVtbtqsCngX8HDg/ybPaxbD2SfLu6Xc4SX0wJKRWVf2EpmvoQpq1g79PEwIH0L3AzRtpFuC5mGZJz9NoBpin/BI4CPgizcp6bwJeVVVfa49/mCYwFtNclezd3r20T/u9n2o/dwbwGGDoi8pIq+OiQ5KkTl5JSJI6GRKSpE6GhCSpkyEhSepkSEiSOhkSkqROhoQkqZMhIUnq9P8BTsJvR2ig+X0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot()\n",
    "\n",
    "sns.countplot(x=kmeans.labels_, ax=ax)\n",
    "ax.set_xlabel('Cluster', fontsize=14)\n",
    "ax.set_ylabel('Galaxy count', fontsize=14)\n",
    "plt.savefig(plot_path + 'cluster_countplot.png', format='png', \n",
    "            bbox_inches='tight', dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting galaxies in each cluster\n",
    "- Loading in raw image data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# loading in pre-processed data\n",
    "data_dir = '../../sdss_CNN/data/cleaned_datasets/'\n",
    "x = np.load(data_dir + 'x.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cl in np.unique(kmeans.labels_):\n",
    "    cluster_evaluate.save_cluster_plot(images = x[kmeans.labels_ == cl], \n",
    "                                       cluster_number = cl,\n",
    "                                       cluster_total= kmeans.n_clusters - 1,\n",
    "                                       plot_path= plot_path) ;\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving cluster labels to file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "label_path = './cluster_tags/'\n",
    "np.save(label_path + 'KMeans-' + str(kmeans.n_clusters) + '_CNN-' + str(cnn_feat) + '_PCA-' + str(pca_feat) + '.npy', \n",
    "        kmeans.labels_)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}