{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "335c108d",
   "metadata": {},
   "source": [
    "# Hands-on clustering \\#1: \n",
    "## Large Multi-spectral image segmentation: MARS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe466679",
   "metadata": {},
   "source": [
    "## 3. Data exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e59c1d62",
   "metadata": {},
   "source": [
    "### 3.1. Read the data and understand the formatting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2939e4a",
   "metadata": {},
   "source": [
    "#### Import first required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c3f78612",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (10, 5)\n",
    "plt.rcParams[\"font.size\"] = 15\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b56fdade",
   "metadata": {},
   "source": [
    "#### 3.1.1 Read the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "350ba4b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "mars = pd.read_csv(\"mars.csv.zip\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5e3237f",
   "metadata": {},
   "source": [
    "#### 3.1.2 Display data summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "af38dd20",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>V1</th>\n",
       "      <th>V2</th>\n",
       "      <th>V3</th>\n",
       "      <th>V4</th>\n",
       "      <th>V5</th>\n",
       "      <th>V6</th>\n",
       "      <th>V7</th>\n",
       "      <th>V8</th>\n",
       "      <th>V9</th>\n",
       "      <th>V10</th>\n",
       "      <th>...</th>\n",
       "      <th>V246</th>\n",
       "      <th>V247</th>\n",
       "      <th>V248</th>\n",
       "      <th>V249</th>\n",
       "      <th>V250</th>\n",
       "      <th>V251</th>\n",
       "      <th>V252</th>\n",
       "      <th>V253</th>\n",
       "      <th>V254</th>\n",
       "      <th>V255</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "      <td>38400.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.373672</td>\n",
       "      <td>0.408012</td>\n",
       "      <td>0.418404</td>\n",
       "      <td>0.420980</td>\n",
       "      <td>0.421814</td>\n",
       "      <td>0.421358</td>\n",
       "      <td>0.426807</td>\n",
       "      <td>0.429181</td>\n",
       "      <td>0.419277</td>\n",
       "      <td>0.418314</td>\n",
       "      <td>...</td>\n",
       "      <td>0.327834</td>\n",
       "      <td>0.358084</td>\n",
       "      <td>0.366237</td>\n",
       "      <td>0.341929</td>\n",
       "      <td>0.376129</td>\n",
       "      <td>0.394448</td>\n",
       "      <td>0.373834</td>\n",
       "      <td>0.370925</td>\n",
       "      <td>0.409842</td>\n",
       "      <td>0.413357</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.098974</td>\n",
       "      <td>0.080131</td>\n",
       "      <td>0.080369</td>\n",
       "      <td>0.081863</td>\n",
       "      <td>0.081315</td>\n",
       "      <td>0.080683</td>\n",
       "      <td>0.080681</td>\n",
       "      <td>0.081414</td>\n",
       "      <td>0.078250</td>\n",
       "      <td>0.077991</td>\n",
       "      <td>...</td>\n",
       "      <td>0.064647</td>\n",
       "      <td>0.056778</td>\n",
       "      <td>0.050779</td>\n",
       "      <td>0.059134</td>\n",
       "      <td>0.052337</td>\n",
       "      <td>0.056993</td>\n",
       "      <td>0.064088</td>\n",
       "      <td>0.062586</td>\n",
       "      <td>0.056333</td>\n",
       "      <td>0.053984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.038469</td>\n",
       "      <td>0.233688</td>\n",
       "      <td>0.259045</td>\n",
       "      <td>0.269580</td>\n",
       "      <td>0.273937</td>\n",
       "      <td>0.274090</td>\n",
       "      <td>0.275866</td>\n",
       "      <td>0.274769</td>\n",
       "      <td>0.273685</td>\n",
       "      <td>0.274631</td>\n",
       "      <td>...</td>\n",
       "      <td>0.143893</td>\n",
       "      <td>0.193508</td>\n",
       "      <td>0.219910</td>\n",
       "      <td>0.174218</td>\n",
       "      <td>0.220138</td>\n",
       "      <td>0.215826</td>\n",
       "      <td>0.198084</td>\n",
       "      <td>0.186781</td>\n",
       "      <td>0.224686</td>\n",
       "      <td>0.219101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>0.299668</td>\n",
       "      <td>0.334763</td>\n",
       "      <td>0.343065</td>\n",
       "      <td>0.345050</td>\n",
       "      <td>0.344962</td>\n",
       "      <td>0.345116</td>\n",
       "      <td>0.350838</td>\n",
       "      <td>0.352490</td>\n",
       "      <td>0.346149</td>\n",
       "      <td>0.345396</td>\n",
       "      <td>...</td>\n",
       "      <td>0.272904</td>\n",
       "      <td>0.312169</td>\n",
       "      <td>0.330923</td>\n",
       "      <td>0.295109</td>\n",
       "      <td>0.339747</td>\n",
       "      <td>0.352853</td>\n",
       "      <td>0.324080</td>\n",
       "      <td>0.321077</td>\n",
       "      <td>0.369166</td>\n",
       "      <td>0.375433</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.365481</td>\n",
       "      <td>0.388543</td>\n",
       "      <td>0.397193</td>\n",
       "      <td>0.399102</td>\n",
       "      <td>0.399245</td>\n",
       "      <td>0.398413</td>\n",
       "      <td>0.403601</td>\n",
       "      <td>0.407299</td>\n",
       "      <td>0.396357</td>\n",
       "      <td>0.394954</td>\n",
       "      <td>...</td>\n",
       "      <td>0.334283</td>\n",
       "      <td>0.362259</td>\n",
       "      <td>0.367263</td>\n",
       "      <td>0.344043</td>\n",
       "      <td>0.375732</td>\n",
       "      <td>0.397163</td>\n",
       "      <td>0.378828</td>\n",
       "      <td>0.372294</td>\n",
       "      <td>0.411261</td>\n",
       "      <td>0.414135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>0.445685</td>\n",
       "      <td>0.487916</td>\n",
       "      <td>0.502658</td>\n",
       "      <td>0.507322</td>\n",
       "      <td>0.508697</td>\n",
       "      <td>0.508175</td>\n",
       "      <td>0.512995</td>\n",
       "      <td>0.515465</td>\n",
       "      <td>0.503711</td>\n",
       "      <td>0.502607</td>\n",
       "      <td>...</td>\n",
       "      <td>0.379558</td>\n",
       "      <td>0.400942</td>\n",
       "      <td>0.398504</td>\n",
       "      <td>0.385884</td>\n",
       "      <td>0.409758</td>\n",
       "      <td>0.433080</td>\n",
       "      <td>0.420058</td>\n",
       "      <td>0.419563</td>\n",
       "      <td>0.447426</td>\n",
       "      <td>0.448774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>0.832065</td>\n",
       "      <td>0.590190</td>\n",
       "      <td>0.580816</td>\n",
       "      <td>0.579686</td>\n",
       "      <td>0.580050</td>\n",
       "      <td>0.579551</td>\n",
       "      <td>0.584240</td>\n",
       "      <td>0.591300</td>\n",
       "      <td>0.569333</td>\n",
       "      <td>0.567649</td>\n",
       "      <td>...</td>\n",
       "      <td>0.559236</td>\n",
       "      <td>0.596519</td>\n",
       "      <td>2.304593</td>\n",
       "      <td>1.393481</td>\n",
       "      <td>0.661885</td>\n",
       "      <td>0.651398</td>\n",
       "      <td>0.624899</td>\n",
       "      <td>0.652349</td>\n",
       "      <td>0.699387</td>\n",
       "      <td>0.655508</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>8 rows × 255 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                 V1            V2            V3            V4            V5  \\\n",
       "count  38400.000000  38400.000000  38400.000000  38400.000000  38400.000000   \n",
       "mean       0.373672      0.408012      0.418404      0.420980      0.421814   \n",
       "std        0.098974      0.080131      0.080369      0.081863      0.081315   \n",
       "min        0.038469      0.233688      0.259045      0.269580      0.273937   \n",
       "25%        0.299668      0.334763      0.343065      0.345050      0.344962   \n",
       "50%        0.365481      0.388543      0.397193      0.399102      0.399245   \n",
       "75%        0.445685      0.487916      0.502658      0.507322      0.508697   \n",
       "max        0.832065      0.590190      0.580816      0.579686      0.580050   \n",
       "\n",
       "                 V6            V7            V8            V9           V10  \\\n",
       "count  38400.000000  38400.000000  38400.000000  38400.000000  38400.000000   \n",
       "mean       0.421358      0.426807      0.429181      0.419277      0.418314   \n",
       "std        0.080683      0.080681      0.081414      0.078250      0.077991   \n",
       "min        0.274090      0.275866      0.274769      0.273685      0.274631   \n",
       "25%        0.345116      0.350838      0.352490      0.346149      0.345396   \n",
       "50%        0.398413      0.403601      0.407299      0.396357      0.394954   \n",
       "75%        0.508175      0.512995      0.515465      0.503711      0.502607   \n",
       "max        0.579551      0.584240      0.591300      0.569333      0.567649   \n",
       "\n",
       "       ...          V246          V247          V248          V249  \\\n",
       "count  ...  38400.000000  38400.000000  38400.000000  38400.000000   \n",
       "mean   ...      0.327834      0.358084      0.366237      0.341929   \n",
       "std    ...      0.064647      0.056778      0.050779      0.059134   \n",
       "min    ...      0.143893      0.193508      0.219910      0.174218   \n",
       "25%    ...      0.272904      0.312169      0.330923      0.295109   \n",
       "50%    ...      0.334283      0.362259      0.367263      0.344043   \n",
       "75%    ...      0.379558      0.400942      0.398504      0.385884   \n",
       "max    ...      0.559236      0.596519      2.304593      1.393481   \n",
       "\n",
       "               V250          V251          V252          V253          V254  \\\n",
       "count  38400.000000  38400.000000  38400.000000  38400.000000  38400.000000   \n",
       "mean       0.376129      0.394448      0.373834      0.370925      0.409842   \n",
       "std        0.052337      0.056993      0.064088      0.062586      0.056333   \n",
       "min        0.220138      0.215826      0.198084      0.186781      0.224686   \n",
       "25%        0.339747      0.352853      0.324080      0.321077      0.369166   \n",
       "50%        0.375732      0.397163      0.378828      0.372294      0.411261   \n",
       "75%        0.409758      0.433080      0.420058      0.419563      0.447426   \n",
       "max        0.661885      0.651398      0.624899      0.652349      0.699387   \n",
       "\n",
       "               V255  \n",
       "count  38400.000000  \n",
       "mean       0.413357  \n",
       "std        0.053984  \n",
       "min        0.219101  \n",
       "25%        0.375433  \n",
       "50%        0.414135  \n",
       "75%        0.448774  \n",
       "max        0.655508  \n",
       "\n",
       "[8 rows x 255 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary = mars.describe()\n",
    "summary\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee6f8866",
   "metadata": {},
   "source": [
    "#### 3.1.3 Data dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c066ab19",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of pixels:   38400\n",
      "Number of channels: 255\n"
     ]
    }
   ],
   "source": [
    "mars_data = mars.to_numpy()\n",
    "print(\"Number of pixels:  \", mars_data.shape[0])\n",
    "print(\"Number of channels:\", mars_data.shape[1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dea62a1",
   "metadata": {},
   "source": [
    "### 3.2. Data preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "925a0b36",
   "metadata": {},
   "source": [
    "#### 3.2.1 Histograms of six randomly selected wave lengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "10f69442",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exploration of different wave length distributions\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22364 characters]",
      "text/plain": [
       "<Figure size 720x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_index = np.random.randint(0, 255, 6)\n",
    "\n",
    "\n",
    "print(\"Exploration of different wave length distributions\")\n",
    "\n",
    "for i in range(6):\n",
    "    #     plt.subplot(2, 3, 1)\n",
    "    #     plt.boxplot(mars_data[:, plot_index])\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    plt.hist(mars_data[:, plot_index[i]], bins=100)\n",
    "    plt.xlabel(\"V%i\" % (plot_index[i] + 1))\n",
    "    plt.ylabel(\"Frequency\")\n",
    "    plt.tight_layout()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f81aaeeb",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "We can see that some wave lengths have clearly different values from different kinds of pixels.\\\n",
    "For some dimensions, they could be used alone to define two distinct clusters easily.\\\n",
    "Overall, the six wave lengths selected appear to be more or less symmetric, although some right spikes are higher than the left spikes, it would not hurt to simply standardize the data.   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da6df892",
   "metadata": {},
   "source": [
    "#### 3.2.2 Box plots of all wave lengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cd8fdaa0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 33064 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boxplot(mars_data)\n",
    "plt.title(\"Box plots of all variables\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25b1af85",
   "metadata": {},
   "source": [
    "##### Interpretations\n",
    "\n",
    "Some dimentions (wave lengths) appear to have much wider spread than others $\\rightarrow$ Need to normalize the data.\\\n",
    "Overall, the median appear to be centered in all the box plots. and they look rather symmetrical $\\rightarrow$ Standardization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "50c866f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 136436 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sample = mars_data[::1000]\n",
    "x = np.arange(0, 255)\n",
    "\n",
    "for i in range(sample.shape[0]):\n",
    "    plt.plot(x, sample[i, :])\n",
    "plt.xlabel(\"Wave length\")\n",
    "plt.ylabel(\"Measured value\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0c68ea6",
   "metadata": {},
   "source": [
    "#### 3.2.4 Data standardization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "aff045da",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "ss = StandardScaler()\n",
    "mars_scaled = ss.fit_transform(mars_data)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7cf6df8d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Visualization of different wave length distributions after standard scaling\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 22388 characters]",
      "text/plain": [
       "<Figure size 720x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Visualization of different wave length distributions after standard scaling\")\n",
    "\n",
    "for i in range(6):\n",
    "    plt.subplot(2, 3, i + 1)\n",
    "    plt.hist(mars_scaled[:, plot_index[i]], bins=100)\n",
    "    plt.xlabel(\"V%i\" % (plot_index[i] + 1))\n",
    "    plt.ylabel(\"Frequency\")\n",
    "    plt.tight_layout()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "765a907f",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "Data are now centered around zero with a spread (standard deviation) of 1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa867370",
   "metadata": {},
   "source": [
    "### 3.3 Dimensionality reduction\n",
    "\n",
    "#### 3.3.1 Conduct PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "58132786",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca = PCA()\n",
    "mars_pca = pca.fit_transform(mars_scaled)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "473e5a13",
   "metadata": {},
   "source": [
    "#### 3.3.2 Visualize explained variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d01d15d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 17220 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8601374108749605\n"
     ]
    }
   ],
   "source": [
    "explained_variance_ratio = pca.explained_variance_ratio_\n",
    "\n",
    "n_bars = 20\n",
    "x = np.arange(1, n_bars + 1)\n",
    "plt.bar(x, explained_variance_ratio[:n_bars])\n",
    "\n",
    "plt.xlabel(\"Number of components\")\n",
    "plt.ylabel(\"Percentage of explained variance\")\n",
    "plt.show()\n",
    "\n",
    "print(sum([explained_variance_ratio[i] for i in range(3)]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afe4acfd",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "The first three PCA components appear to be explaining most of the variance in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15ef4729",
   "metadata": {},
   "source": [
    "#### 3.3.3 Variable factor map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4290cb31",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 152900 characters]",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "coord1 = pca.components_[0] * np.sqrt(pca.explained_variance_[0])\n",
    "coord2 = pca.components_[1] * np.sqrt(pca.explained_variance_[1])\n",
    "\n",
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "for i, j, nom in zip(coord1, coord2, mars.columns):\n",
    "    plt.text(i, j, nom, fontsize=10)\n",
    "    plt.arrow(0, 0, i, j, color=\"r\", width=0.0001)\n",
    "\n",
    "plt.axis((-1, 1, -1, 1))\n",
    "plt.gcf().gca().add_artist(plt.Circle((0, 0), radius=1, color=\"b\", fill=False))\n",
    "\n",
    "plt.title(\"Variables factor map - PCA\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3119be95",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "Most of the wave lengths have projection on the first two PCA dimensions with norm 1.\\\n",
    "Most of them are represented entirely by component 1, other entirely by component 2.\\\n",
    "Some isolated wave lengths are not well represented in the first two PCA dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26955b8b",
   "metadata": {},
   "source": [
    "#### 3.3.4 Individual factor maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9005489b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 190516 characters]",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(10, 10))\n",
    "\n",
    "plt.scatter(mars_pca[:, 0], mars_pca[:, 1], s=1, linewidths=1)\n",
    "plt.title(\"Individuals factor map - PCA\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efe72fd5",
   "metadata": {},
   "source": [
    "##### Interpretations\n",
    "\n",
    "It is hard to distinguish visually the five groups identified by experts.\\\n",
    "Two, three groups appear to be clearly present in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "914979de",
   "metadata": {},
   "source": [
    "#### 3.3.5 Reduce dimension to three"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f04b7c35",
   "metadata": {},
   "outputs": [],
   "source": [
    "mars_reduced = mars_pca[:, :3]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25c25c3a",
   "metadata": {},
   "source": [
    "## 4. Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1b8020d",
   "metadata": {},
   "source": [
    "### 4.1. K-means"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0499c35a",
   "metadata": {},
   "source": [
    "#### Additional imports "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d1fe2f10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans, AgglomerativeClustering\n",
    "from sklearn.metrics import silhouette_score\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "650aba6a",
   "metadata": {},
   "source": [
    "#### 4.1.1 Choose the number of clusters\n",
    "\n",
    "Either use the elbow method or the silhouette method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edbac1f5",
   "metadata": {},
   "source": [
    "#### 4.1.2 Elbow method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1df27f92",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12900 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "inertia = []\n",
    "for k in range(1, 11):\n",
    "    kmeans = KMeans(n_clusters=k)\n",
    "    kmeans.fit(mars_reduced)\n",
    "    inertia.append(kmeans.inertia_)\n",
    "inertia = np.array(inertia)\n",
    "\n",
    "plt.scatter(range(2, 11), inertia[1:])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1470a547",
   "metadata": {},
   "source": [
    "#### Interpretation\n",
    "\n",
    "Choose K=5"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaec9ce2",
   "metadata": {},
   "source": [
    "#### 4.1.3 Apply K-means on PCA reduced data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "b2e92da5",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans_pca = KMeans(n_clusters=5, random_state=0)\n",
    "clusters_pca = kmeans_pca.fit_predict(mars_reduced)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a825e1f",
   "metadata": {},
   "source": [
    "#### 4.1.4 Display results on individual factor map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6500d1f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 173856 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(mars_pca[:, 0], mars_pca[:, 1], c=clusters_pca, s=1, linewidths=1)\n",
    "plt.title(\"Individuals factor map - PCA\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73e7aa0d",
   "metadata": {},
   "source": [
    "#### 4.1.5 Display the image of MARS geological composition found by clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "12093013",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30896 characters]",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_pixel_x = 300\n",
    "n_pixel_y = 128\n",
    "\n",
    "mars_image = clusters_pca.reshape((n_pixel_x, n_pixel_y))\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.imshow(mars_image, interpolation=\"nearest\", aspect=\"auto\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42036349",
   "metadata": {},
   "source": [
    "#### 4.1.6 Curve representing the wave lengths of the centroids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f1ceca9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 88068 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(0, 255)\n",
    "for i in range(5):\n",
    "    plt.plot(x, np.mean(mars_data[clusters_pca == i], axis=0))\n",
    "plt.xlabel(\"Wave length\")\n",
    "plt.ylabel(\"Measured value\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a6f785c",
   "metadata": {},
   "source": [
    "#### 4.1.7 Repeat the process with complete data and compare the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "9f9ea693",
   "metadata": {},
   "outputs": [],
   "source": [
    "kmeans = KMeans(n_clusters=5, random_state=0)\n",
    "clusters_full = kmeans.fit_predict(mars_data, 5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf958021",
   "metadata": {},
   "source": [
    "#### Import and use external metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0a910a4d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  49,  662,  118,    0, 4331],\n",
       "       [1699,    0,    0, 8474,    0],\n",
       "       [   0, 8361,  866,    0,   50],\n",
       "       [4527,    0,    8,  748,  958],\n",
       "       [   0,   66, 7483,    0,    0]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "confusion_matrix(clusters_full, clusters_pca)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "8075c617",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Normalized Mutual Information score: 0.7420448071402478\n",
      "Fowlkes Mallows score:  0.7830336140166589\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import normalized_mutual_info_score, fowlkes_mallows_score\n",
    "\n",
    "similarity_nmi = normalized_mutual_info_score(clusters_full, clusters_pca)\n",
    "similarity_fm = fowlkes_mallows_score(clusters_full, clusters_pca)\n",
    "\n",
    "print(\"Normalized Mutual Information score:\", similarity_nmi)\n",
    "print(\"Fowlkes Mallows score: \", similarity_fm)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7617905b",
   "metadata": {},
   "source": [
    "##### Interpretations\n",
    "\n",
    "The scores obtained with both external metrics are high, meaning that the classifications obtained are similar with and without PCA."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ec5fd6c",
   "metadata": {},
   "source": [
    "### 4.2 Agglomerative Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6871c039",
   "metadata": {},
   "source": [
    "#### 4.2.1. Choose the number of clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b08d8309",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 10640 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Take only a fraction of the initial data to debug.\n",
    "mars_reduced_samples = mars_reduced[::10]\n",
    "\n",
    "ac = AgglomerativeClustering(n_clusters=5, linkage=\"ward\", compute_distances=True)\n",
    "clusters_ac = ac.fit_predict(mars_reduced_samples)\n",
    "\n",
    "distances = ac.distances_\n",
    "\n",
    "n_sizes = 20\n",
    "x = np.arange(n_sizes, 0, -1)\n",
    "y = ac.distances_[-20:]\n",
    "\n",
    "plt.scatter(x, y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "108f9bd5",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "In this plot, the first point is much higher than the other, making the y-scale of all the interesting points too small and making it hard to discover the elbow. $\\rightarrow$ Remove the first point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "866ae33c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 11148 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_sizes = 10\n",
    "x = np.arange(n_sizes, 1, -1)\n",
    "y = ac.distances_[-n_sizes:-1]\n",
    "\n",
    "plt.scatter(x, y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7c28be7",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "In these plots, either three or five appear to be a smart choice for the number of clusters. It does not allow us to decide precisely what is the right number (maybe it does not even exist). However, it gives us the right range compared to what experts said (not 20, not 100). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53314ab4",
   "metadata": {},
   "source": [
    "#### 4.2.2. Visualize agglomerative clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "72dd0f29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fa7ce44acd0>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 113296 characters]",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(\n",
    "    mars_reduced_samples[:, 0],\n",
    "    mars_reduced_samples[:, 1],\n",
    "    c=clusters_ac,\n",
    "    s=1,\n",
    "    linewidths=1,\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edcd5565",
   "metadata": {},
   "source": [
    "## 5. Comparison of clustering algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c032c99",
   "metadata": {},
   "source": [
    "### 5.1. Compare different approaches"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26eb1227",
   "metadata": {},
   "source": [
    "#### Load experts ground truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "86f61b7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters_expert = pd.read_csv(\"mask.csv\")\n",
    "clusters_expert = clusters_expert.to_numpy()\n",
    "clusters_expert = np.reshape(clusters_expert, 38400)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd42f945",
   "metadata": {},
   "source": [
    "#### Compare different approaches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "056ad896",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results for K-means applied on first three dimensions of PCA\n",
      "FM : 0.6015665557841873\n",
      "NMI: 0.5387823662410469\n",
      "Results for K-means applied on original data (normalized)\n",
      "FM : 0.6263019384432416\n",
      "NMI: 0.5582715747398275\n",
      "Results for Agglomerative Clustering applied on original data (small subset, normalized)\n",
      "FM : 0.6276636966595519\n",
      "NMI: 0.5378671553934923\n"
     ]
    }
   ],
   "source": [
    "print(\"Results for K-means applied on first three dimensions of PCA\")\n",
    "print(\"FM :\", fowlkes_mallows_score(clusters_pca, clusters_expert))\n",
    "print(\"NMI:\", normalized_mutual_info_score(clusters_pca, clusters_expert))\n",
    "\n",
    "print(\"Results for K-means applied on original data (normalized)\")\n",
    "print(\"FM :\", fowlkes_mallows_score(clusters_full, clusters_expert))\n",
    "print(\"NMI:\", normalized_mutual_info_score(clusters_full, clusters_expert))\n",
    "\n",
    "clusters_expert_samples = clusters_expert[::10]\n",
    "\n",
    "print(\n",
    "    \"Results for Agglomerative Clustering applied on original data (small subset, normalized)\"\n",
    ")\n",
    "print(\"FM :\", fowlkes_mallows_score(clusters_ac, clusters_expert_samples))\n",
    "print(\"NMI:\", normalized_mutual_info_score(clusters_ac, clusters_expert_samples))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6080dc4b",
   "metadata": {},
   "source": [
    "#### Interpretation\n",
    "\n",
    "Results appear relatively similar.\\\n",
    "There is some randomness in both versions of the algorithms (Random initialization for K-means, random subset of data for agglomerative clustering).\\\n",
    "To obtain more relevant results, several runs should be conducted and averaged for K-means and all the data should be used for agglomerative clustering."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e168c94e",
   "metadata": {},
   "source": [
    "#### Several runs for K-means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "58ae2714",
   "metadata": {},
   "outputs": [],
   "source": [
    "FMs = []\n",
    "NMIs = []\n",
    "for i in range(10):\n",
    "    clusters = KMeans(5).fit_predict(mars_reduced)\n",
    "    FMs.append(fowlkes_mallows_score(clusters, clusters_expert))\n",
    "    NMIs.append(normalized_mutual_info_score(clusters, clusters_expert))\n",
    "FMs = np.array(FMs)\n",
    "NMIs = np.array(NMIs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "f27cc08c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results for K-means applied on first three dimensions of PCA\n",
      "FM : 0.602208 ± 0.000842\n",
      "NMI: 0.539388 ± 0.001018\n"
     ]
    }
   ],
   "source": [
    "print(\"Results for K-means applied on first three dimensions of PCA\")\n",
    "print(\"FM : %f \\u00B1 %f\" % (FMs.mean(), FMs.std()))\n",
    "print(\"NMI: %f \\u00B1 %f\" % (NMIs.mean(), NMIs.std()))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9b39c5d",
   "metadata": {},
   "source": [
    "##### Interpretation\n",
    "\n",
    "Not so much variability in the end."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c00b970",
   "metadata": {},
   "source": [
    "### 5.2 Try Gaussian Mixture Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "5bf0c2b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Results for K-means applied on original data (normalized)\n",
      "FM : 0.5926714304424104\n",
      "NMI: 0.5585575432531626\n"
     ]
    }
   ],
   "source": [
    "from sklearn.mixture import GaussianMixture\n",
    "\n",
    "clusters_gmm = GaussianMixture(5).fit_predict(mars_reduced)\n",
    "\n",
    "print(\"Results for K-means applied on original data (normalized)\")\n",
    "print(\"FM :\", fowlkes_mallows_score(clusters_gmm, clusters_expert))\n",
    "print(\"NMI:\", normalized_mutual_info_score(clusters_gmm, clusters_expert))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b44b7d4",
   "metadata": {},
   "source": [
    "##### Interpretations\n",
    "\n",
    "Not better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92987c94",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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
}
