GMM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm, multivariate_normal\n",
    "import scipy.io as spio\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.metrics import adjusted_rand_score, adjusted_mutual_info_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read-in, visualise, convert, flatten and concatenate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "IP_data_mat=spio.loadmat(\"data/IP/Indian_pines_corrected.mat\")\n",
    "IP_gt_mat=spio.loadmat(\"data/IP/Indian_pines_gt.mat\")\n",
    "IP_image=IP_data_mat[\"indian_pines_corrected\"]\n",
    "IP_gt=IP_gt_mat['indian_pines_gt']\n",
    "\n",
    "pavia_mat = spio.loadmat('data/PaviaU/PaviaU.mat')\n",
    "pavia_gt_mat = spio.loadmat('data/PaviaU/PaviaU_gt.mat')\n",
    "pavia_image = pavia_mat['paviaU']\n",
    "pavia_gt = pavia_gt_mat['paviaU_gt']\n",
    "\n",
    "salinas_mat = spio.loadmat('data/Salinas/salinas.mat')\n",
    "salinas_gt_mat = spio.loadmat('data/Salinas/salinas_gt.mat')\n",
    "salinas_image = salinas_mat['salinasA_corrected']\n",
    "salinas_gt = salinas_gt_mat['salinasA_gt']\n",
    "\n",
    "n, m, c = salinas_image.shape\n",
    "salinas_image = np.reshape(salinas_image, (-1, c))\n",
    "salinas_gt_flat = np.reshape(salinas_gt, -1)\n",
    "salinas_image = salinas_image[salinas_gt_flat!=0]\n",
    "\n",
    "# PCA\n",
    "n_features = 5\n",
    "pca = PCA(n_components=n_features)\n",
    "pCs_salinas = pca.fit_transform(salinas_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GMM:\n",
    "    \"\"\"\n",
    "    Class implementing a Gaussian Mixture Model (GMM) with multivariate gaussians\n",
    "        Requires: potential_features (data), n_coms (number of clusters), em_steps(number of iterations of EM) and\n",
    "        n_features (the number of features to extract from potential_features)\n",
    "        \n",
    "        Contains: \n",
    "            GMM.run(), run EM with parameters specified at initialisation, will create appropriate plots\n",
    "    \"\"\"\n",
    "    def __init__(self, potential_features, n_coms=2, em_steps=10, n_features=2):\n",
    "        # Select desired features\n",
    "        self.X = potential_features[:, :n_features]\n",
    "        # Number of steps em will run for\n",
    "        self.em_steps = em_steps\n",
    "        # Number of components/clusters\n",
    "        self.n_coms = n_coms\n",
    "        self.mu = None\n",
    "        self.pi = None\n",
    "        self.cov = None\n",
    "        self.XY = None\n",
    "        self.n_features=n_features\n",
    "    def run(self):\n",
    "        # Regulariser for the covariance matrix\n",
    "        self.reg_cov = 1e-4*np.identity(len(self.X[0]))\n",
    "        \n",
    "        #x,y = np.meshgrid(np.sort(self.X[:,0]),np.sort(self.X[:,1]))\n",
    "        #self.XY = np.array([x.flatten(),y.flatten()]).T\n",
    "        \n",
    "        # Initialise mean value size (Number of components, number of features) and make float\n",
    "        chosen = np.random.choice(len(self.X), self.n_coms, replace = False)\n",
    "        self.mu = self.X[chosen]*1.0\n",
    "        \n",
    "        # Initialise zero matrices of size (number of components, number of features, number of features)\n",
    "        shape = self.n_coms, len(self.X[0]), len(self.X[0])\n",
    "        self.cov = np.full(shape, np.cov(self.X, rowvar = False))\n",
    "        #self.cov = np.sqrt(self.cov) THIS WAS BREAKING THE COIN PIXELS\n",
    "        self.cov += self.reg_cov\n",
    "        \n",
    "        initial_mu = self.mu.copy()\n",
    "        initial_cov = self.cov.copy()\n",
    "        \n",
    "        # Initialise uniform pi values\n",
    "        self.pi = np.full(self.n_coms, 1/self.n_coms)\n",
    "        \n",
    "        log_likelihoods = []\n",
    "        \n",
    "        # EM Loop\n",
    "        i=0\n",
    "        while True:\n",
    "            # E Step\n",
    "            #print(\"EM step %d/%d.\"%(i+1, self.em_steps), end='\\r')\n",
    "            \n",
    "            # Compute responsibilities\n",
    "            r_ic = np.zeros((len(self.X), self.n_coms))\n",
    "            for c in range(self.n_coms):\n",
    "                prior = self.pi[c]\n",
    "                #print(self.mu[c], self.cov[c])\n",
    "                likelihood = multivariate_normal(self.mu[c], self.cov[c]).pdf(self.X)\n",
    "                r_ic[:, c] = prior * likelihood\n",
    "            r_ic = r_ic/(r_ic.sum(axis = 1, keepdims = 1)+1e-9)\n",
    "        \n",
    "            # M Step\n",
    "            \n",
    "            self.pi = []\n",
    "            self.mu = []\n",
    "            self.cov = []\n",
    "            \n",
    "            for c in range(self.n_coms):\n",
    "                # Update mean of each component\n",
    "                m_c = np.sum(r_ic[:, c], axis=0)\n",
    "                mu_c = (1/m_c)*np.sum(self.X*r_ic[:,c].reshape(len(self.X), 1), axis=0)\n",
    "                self.mu.append(mu_c)\n",
    "                \n",
    "                # Update covariance matrix of each component\n",
    "                cov_c = (1/m_c)*np.dot((np.array(r_ic[:,c]).reshape(len(self.X),1)*(self.X-mu_c)).T,(self.X-mu_c))\n",
    "                cov_c += self.reg_cov\n",
    "                self.cov.append(cov_c)\n",
    "                \n",
    "                # Update pi values of each component\n",
    "                self.pi.append(m_c/np.sum(r_ic))\n",
    "            \n",
    "            self.mu = np.array(self.mu)\n",
    "            self.cov = np.array(self.cov)\n",
    "            # Save log likelihood of each iteration\n",
    "            log_likelihoods.append(np.log(np.sum([k*multivariate_normal(self.mu[i],self.cov[j]).pdf(self.X) for k,i,j in zip(self.pi,range(len(self.mu)),range(len(self.cov)))])))\n",
    "            i+=1\n",
    "            print(np.round(log_likelihoods[i-1], 2), i, end='\\r')\n",
    "            if i>1:\n",
    "                if np.abs(log_likelihoods[i-2]-log_likelihoods[i-1])<1e-4:\n",
    "                    break\n",
    "\n",
    "        print('Training complete')\n",
    "        plt.plot(range(0, i, 1),log_likelihoods)\n",
    "        plt.xlabel('Iterative steps')\n",
    "        plt.ylabel('Log-likelihood')\n",
    "        plt.savefig('log_likelihood_graph.png')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training complete\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8XXX9x/HXJ6Nt6Eg3dNIWWqBAU0qggKBYAbEoW/byR6mgqKA/ZIgK/H4qAor158AyRcqQUWUvAZHRQtKRlg5aOtN0pCNpOtJmfH5/3BO4DTfJTXJvTm7u+/l43Mc943vO+Ryu5tPvON9j7o6IiEhrZYQdgIiIdAxKKCIikhBKKCIikhBKKCIikhBKKCIikhBKKCIikhBKKCIikhBKKCIikhBKKCIikhBZYQfQlvr27evDhg0LOwwRkZRSWFi40d37NVUurRLKsGHDKCgoCDsMEZGUYmYr4ymnJi8REUkIJRQREUkIJRQREUkIJRQREUkIJRQREUkIJRQREUkIJRQREUmItHoORSQVVdXUUr6zivKdVVRUVrO7upbqmlqqap2a2lqqapzqGqc6WK7b5u7U1Dq1DrXuePBdt14b7HMirwGvexv4Hi8Fj/WKcDPss0Xq1iLLERkZtsf+un3R6wAZFr0vWA7O39ixe5Spd87Gjoc9rxfZ8tl5ghIQo0ys+/10X/SxUWX3PDbGcfXKEuP4+mWjzxV7X/Rxn63s23cvenTJJpmUUERCtmFrJQ+9t4KN23Z9mjjKdlSxNVjevrsm1Pii/0DFyi+SGh761hEcf0D/pF4jlIRiZk8ABwSrPYEydx9br8wBwBNRm0YAP3P335nZLcAVQGmw7yZ3fzG5UYsk3q7qGq54uID5JVvp07UTPffKJjcnm8G99iJ3YGS5bltuTjY9crLolJlJVqaRlWFkZWaQlWFkZ2aQlWlkZ2R8ui8jw8gwI9MMy4jUBjLss1pBphkWbKsT/S/fpnhQ64FIrcaDlbpaz6f7otbrynmwnbraU73teyzXP9abOOce26PK1YsrVmx87lx1d0e9MnuWZY/yn9X66n1F/ff6rHCsmqFHnexzOTzWdWLFUO/AQwbl1j9TwoWSUNz93LplM/sNUB6jzGJgbFAmE1gDTI8qcre735XkUEWS6n+eX8Dc4nLuuehwTj5kn7DDaZa6JqqoLWGFIu1EqE1eFvnn0DnAhCaKfgX4xN3jmk9GJBVMn13MIzNW8e0vjki5ZCISS9ijvI4D1rv7kibKnQc8Vm/b1WZWZGYPmFmv5IQnkhyL1m3lxmfmceTw3lz31QOaPkAkBSQtoZjZ62Y2P8bntKhi5/P5RFH/PJ2AU4Enozb/GdiPSJPYWuA3jRw/2cwKzKygtLS0oWIibaaisoqrHplF9y7Z/OGCw8jKDPvfdSKJkbQmL3c/obH9ZpYFnAkc3sSpvgbMcvf1Uef+dNnM7gWebySOqcBUgPz8fI1RkVC5O9c9WcSqzTt4dNJ4+nfvEnZIIgkT5j+NTgAWuXtxE+U+V4sxswFRq2cA8xMcm0hS3P/Ocl7+aB03nHwg40f0CTsckYQKM6F8rl/EzAaa2YtR63sBJwLP1Dv2DjObZ2ZFwJeBa5MdrEhrfbB8M796aREnH7wPk44bHnY4IgkX2igvd78sxrYSYGLU+g7gc/+Mc/eLkxqcSIJtqKjku4/OYmjvvbjzm2Oa9byHSKrQk/IiSVZdU8vVj86morKKv11+JN2TPP2FSFiUUESS7M5XF/PB8s3cfW4eB+7TI+xwRJJG4xVFkujl+ev4y7+XcdFRQznjsMFhhyOSVEooIkny8foKrntyLnmDc/np10eHHY5I0qnJSyTBNm/fzf+9sYRHZqyke5ds/njhODpnZYYdlkjSKaGIJEhlVQ0PvruCP721lO27qjknfwjXnjiKvXvo4UVJD0ooIq1UW+v8Y84a7nplMSXllUw4sD83fO1ARu3dPezQRNqUEopIK7yzZCO/fHEhC9Zu5ZBBPbjrnDyO2a9v2GGJhEIJRaQFFq3byq9eXMS/Py5lUM8cppw3lm+MGfjpq29F0pESikgz7Nxdw+0vLeRvM1bSrXMWN008kEuOHkaXbHW6iyihiMTpo5JyfvD4HJZu2MZlxwzjmhNG0nOvTmGHJdJuKKGINKG21nng3eXc8fJieu6VzSOXj+fYkeonEalPCUWkERu2VvKjJ+fynyUbOXH03vz6rDH07qpaiUgsSigiDXhtwXp+/NRcdlbV8MszDuX8I4dolmCRRiihiNSzc3cN//vCAqbNXMXBA3sw5bzD2L9/t7DDEmn3lFBEokR3vE/+4gh+dNIoTZsiEiclFBEi73p/8N0V3P7SInW8i7SQEoqkvcqqGm58Zh7TZ6/hhIP25o6z1fEu0hJKKJLW1pbv5Nt/K6SouJwfnTiKqyfsr453kRZSQpG0VbhyM9/+2ywqq2q495J8Thy9d9ghiaS00F6wZWZjzWyGmc0xswIzO7KBcpea2ZLgc2nU9sPNbJ6ZLTWz35v+WSnN8MSHqzhv6gy6dc5k+neOUTIRSYAw39h4B3Cru48Ffhas78HMegM/B8YDRwI/N7Newe4/A5OBkcHn5LYIWlJbVU0tP//nfK5/eh5HjejDP797LCM1zbxIQoTZ5OVAj2A5FyiJUearwGvuvhnAzF4DTjazt4Ae7v5+sP1h4HTgpWQHLalr8/bdfGdaITOWbWbyF0fw468eQFam3oItkihhJpRrgFfM7C4iNaVjYpQZBKyOWi8Otg0Klutv/xwzm0ykJsPQoUNbH7WkpAUlW7ni4QJKt+3i7nPzOOOwwWGHJNLhJDWhmNnrwD4xdv0E+Apwrbs/bWbnAPcDJ9Q/RYxjvZHtn9/oPhWYCpCfnx+zjHRsz80t4cdPFZGbk82T3z6avCE9ww5JpENKakJx9/oJ4lNBM9UPgtUngftiFCsGjo9aHwy8FWwfXG97rCYzSWNLN2zjf19YwFuLSxk3tCf3XHw4/bvr/e4iyRJmk1cJ8CUiCWICsCRGmVeAX0Z1xJ8E3Ojum82swsyOAmYClwD/l/yQJRWU76hiyr+W8PD7K8jJzuQnEw/i0mOG0SlL/SUiyRRmQrkCmGJmWUAlQT+HmeUDV7r7pCBx/A/wYXDMbXUd9MBVwENADpHOeHXIp7maWuexD1bxm1cXU7azivOOGMqPThpF326dww5NJC2Ye/p0K+Tn53tBQUHYYUgSvPfJRm57bgGL1lUwfnhvfvaN0Rw8MDfssEQ6BDMrdPf8psrpSXlJaas27eCXLy7k5Y/WMbhXDn+6cBxfO2QfTZ8iEgIlFElJO3fX8Ic3l3Dvf5aTlWH890mjmHTcCLpka6p5kbAooUjKeWfJRm6aPo9Vm3dw5mGD+PHJB7JPrkZviYRNCUVSxpbtu/nFiwt5qrCY4X278vjkozhqRJ+wwxKRgBKKtHvuznNFa7n12Y8o31nFd7+8H9+bMFLNWyLtjBKKtGtrynby03/M541FG8gbnMsjk8Zz0IAeTR8oIm1OCUXapZpa52/vr+DOVxZT6/DTr4/msmOGkZmh0Vsi7ZUSirQ7i9dVcMMzRcxeVcYXR/XjF6cfwpDee4Udlog0QQlF2pUZyzZx8f0z6d4lm9+dO5bTxg7UMyUiKUIJRdqNnbtruP7pIgb2zOGZq46hj6ZMEUkpSijSbvzu9Y9ZuWkHj14xXslEJAVp+lVpF4qKy7j3P8s4/8ghHLNf37DDEZEWUEKR0FXV1PLjp4ro260zN3ztoLDDEZEWUpOXhG7q28tYtK6CqRcfTm5OdtjhiEgLqYYioVq6YRtTXl/CKYcO4KSDY70tWkRShRKKhKa21rnxmSJyOmVyy6kHhx2OiLSSEoqEZtrMlXy4Ygs//fpo+nXXqC6RVKeEIqFYU7aT219axHEj+3LWuEFhhyMiCaCEIm3O3bl5+jxqHX55xqF6El6kgwgloZjZWDObYWZzzKzAzI5soMz7ZvaRmRWZ2blR+x4ys+XB8XPMbGzb3oG0xrNzS3hzcSnXffUAzdEl0oGENWz4DuBWd3/JzCYG68fXK7MDuMTdl5jZQKDQzF5x97Jg/3Xu/lTbhSyJsGnbLm559iPGDunJpccMCzscEUmgsBKKA3UvtcgFSj5XwP3jqOUSM9sA9APK6peV1HHb8wvYtquaO84eo6noRTqYsPpQrgHuNLPVwF3AjY0VDprEOgGfRG3+RdAUdreZaYhQCnhj0Xr+OaeE7355f0bt3T3scEQkwZKWUMzsdTObH+NzGnAVcK27DwGuBe5v5DwDgL8B33L32mDzjcCBwBFAb+D6Ro6fHPTTFJSWlibo7qS5KiqruHn6fEbt3Y3vHL9/2OGISBKYu7f9Rc3KgZ7u7hYZ4lPu7p97r6uZ9QDeAn7l7k82cK7jgf929683dd38/HwvKChoVezSfNt3VTPprwXMXL6Jp646hnFDe4Udkog0g5kVunt+U+XCavIqAb4ULE8AltQvYGadgOnAw/WTSVBrIUhGpwPzkxqttNjWyioueeADZi7fxG/OyVMyEenAwuqUvwKYYmZZQCUwGcDM8oEr3X0ScA7wRaCPmV0WHHeZu88BpplZP8CAOcCVbRy/xGHL9t1c8sAHLFy7lT9cMI6Jhw4IOyQRSaJQmrzCoiavtlNasYuL75/Jso3bueeicUw4cO+wQxKRFoq3yavRGoqZPUdkiG9M7n5qC2KTDm5t+U4uvG8ma8sqeeDSIzh2pF6YJZIOmmryuiv4PhPYB3gkWD8fWJGkmCSFrd68gwvum8GW7VU8fPmRHDGsd9ghiUgbaTShuPu/Aczsf9z9i1G7njOzt5MamaScZaXbuPC+mezYXcO0SePJG9Iz7JBEpA3FO8qrn5mNqFsxs+FEnloXAWDxugrO+csMdlfX8tgVRymZiKSheEd5XQu8ZWbLgvVhBCOzROavKefi+2fSKSuDaZOOYv/+egpeJB3FlVDc/WUzG0nk6XSARe6+K3lhSSpwd95aXMr3H59Njy7ZPHrFePbt0zXssEQkJHElFDPLBr5N5LkQiNRW/uLuVUmLTNq1uavL+PXLi3jvk03s168rD18+nkE9c8IOS0RCFG+T15+BbOBPwfrFwbZJyQhK2q/lG7dz1yuLeWHeWnp37cTPvzGaC8YPpXNWZtihiUjI4k0oR7h7XtT6G2Y2NxkBSfu0YWslU/61hMc/XE3nrAy+/5WRXHHccLp3yQ47NBFpJ+JNKDVmtp+7fwIQjPiqSV5Y0l5UVFYx9e1l3Pef5VTV1HLh+KF8b8JI+nXXGwNEZE/xJpTrgDeDUV4G7At8K2lRSeh2VdfwyIxV/OGNJWzZUcU38gbyoxNHMayvOt1FJLZ4R3n9KxjldQCRhKJRXh3U0g0VPPHhap6ZtYZN23dz3Mi+XH/ygRwyKDfs0ESkndMoL2HH7mqeL1rL3z9cTcHKLWRlGCcctDeXHL0vx+yvebhEJD4a5ZWm3J25xeU88eFqnptbwrZd1ezXrys3TTyQM8cNpm839ZGISPNolFea2bJ9N/+Ys4YnPlzNonUV5GRncsqYAZx3xBAO37cXkXeWiYg0n0Z5pYHyHVW8smAdLxSt5d2lG6mudfIG5/LLMw7lG3kDNPRXRBJCo7w6qPIdVby6YB0vzFvLO0siSWRI7xwuP244p+UNYvTAHmGHKCIdjEZ5dSDlO6t4bcF6Xigq4Z2lG6mqcQb3yuHyY4dzypgBHDooV01aIpI0zXmn/OFEZhnOAvLMDHd/OClRSdzcnfc+2cRD763grcUbqKpxBvXM4VtfGM4phw5gzGAlERFpG/EOG/4bsB8wh8/6ThxQQgnJruoanp1Twv3vLGfRugr6duvEpUcP4+t5A8lTEhGREMRbQ8kHRrt7g++Xby4zGwvcA3QBqoHvuPsHMcrVAPOC1VV177EPXvL1ONAbmAVc7O67ExVfe7Vp2y6mzVzFw++vZOO2XRywd3fuOHsMp+YNpEu2JmgUkfDEm1DmE3mn/NoEXvsO4FZ3f8nMJgbrx8cot9Pdx8bY/mvgbnd/3MzuAS4n8mxMh7RkfQUPvLucZ2atYVd1Lccf0I9Jx47gC/v3UW1ERNqFRhOKmT1HpGmrO7DAzD4APu2Mr6sttJADdUONcoGSeA+0yF/QCcAFwaa/ArfQARPKe0s38pe3l/Hvj0vpnJXBmeMGc/mxw/RWRBFpd5qqodyVxGtfA7xiZncRebf9MQ2U62JmBUSaxW53938AfYAyd68OyhQDg5IYayienVvC9x+bTb/unfnRiaO48Kh96d21U9hhiYjE1GhCcfd/t+bkZvY6kaay+n4CfAW41t2fNrNzgPuBE2KUHeruJcHDlG+Y2Txga6xwG4hhMjAZYOjQoS24i3As3VDBDU8Xcfi+vZg2abz6R0Sk3bPG+tnN7B13P9bMKtjzD7YB7u4tfjrOzMqBnu7uQRNWeVPnM7OHgOeBp4FSYB93rzazo4Fb3P2rjR2fn5/vBQUFLQ25zezYXc1pf3iXzdt388L3j2Of3C5hhyQiaczMCt09v6lyGY3tdPdjg+/u7t4j6tO9NckkUAJ8KVieACypX8DMeplZ52C5L/AFYEEw2uxN4Oyg6KXAP1sZT7vg7tz0zDyWlm5jynmHKZmISMpoqlO+d2P73X1zK659BTDFzLKASoJmKTPLB65090nAQcBfzKyWSPK73d0XBMdfDzxuZv8LzCbSZJbyHv1gFf+YU8IPTxzFsSM1dbyIpI6mmryWE2nqijUu1d19RLICS4b23uQ1r7ics/78Hkfv14cHLzuCjAwNBxaR8MXb5NVUp/zwxIUkjSnfUcVV0wrp260Tvzt3rJKJiKScRvtQ6ljERWb202B9qJkdmdzQ0kdtrfPDv89h/dZK/njhOHppaLCIpKC4EgqRNzUezWcPElYAf0xKRGnonrc/4V+LNnDzKaM5bGivsMMREWmReKdeGe/u48xsNoC7bzEz/TM6Ad7/ZBN3vbKYr48ZwCVH7xt2OCIiLRZvDaXKzDIJnkUxs35AbdKiShMbtlbyvcdmM6xvV24/a4zm5BKRlBZvQvk9MB3ob2a/AN4BfpW0qNJAdU0tVz82m+27qrnnosPp1rk5r6YREWl/4n1j4zQzKyQyXYoBp7v7wqRG1sHd9erHfLB8M789J49Re2uiRxFJffG+YOtyd78fWBS17XZ3vyFpkXVgc1eXcc+/P+GC8UM5c9zgsMMREUmIeNtZzjazSnefBmBmfwI6Jy+sju0/S0oB+PFXDwg5EhGRxIk3oZwJPBtMgfI1YLO7fyd5YXVshSu3MLJ/N3rupYFyItJxNNopb2a9g/m8coBJwI+JTB1/W1PzfElstbXOrFVljNPzJiLSwTRVQynks7m86r5PCT4OpNRcXu3Bso3bKd9ZxeH7KqGISMeiubza2KyVWwAYt2/PkCMREUmspqavn+Dub5jZmbH2u/szyQmr45q1agu5OdmM6Nst7FBERBKqqSavLwFvAN+Isc8BJZRmKly5hXFDe2o2YRHpcJpq8vp58P2ttgmnYyvfUcWSDds4NW9g2KGIiCRcU01eP2xsv7v/NrHhdGyzV0f6T9QhLyIdUVNNXpoTJIFmrdxChkHeEHXIi0jH01ST161tFUg6KFy1hQP36UFXTQQpIh1QvLMNf8rMZiUjkI6uptaZs6pMzV0i0mE1O6EQebixVcxsrJnNMLM5ZlYQ63XCZvblYH/dp9LMTg/2PWRmy6P2jW1tTMm2eF0F23fXKKGISIfVkraXFxJw3TuAW939JTObGKwfH13A3d8ExkJkChhgKfBqVJHr3P2pBMTSJmatCh5o1JQrItJBNbuG4u43J+C6DvQIlnOBkibKnw285O47EnDtUMxauYW+3TozpHdO2KGIiCRFXAnFzCrMbGu9z2ozm25mLZnP6xrgTjNbDdwF3NhE+fOAx+pt+4WZFZnZ3WbW4FT6ZjY5aFYrKC0tbUGoiVG4KvJAo17zKyIdVbw1lN8C1wGDgMHAfwP3Ao8DD8Q6wMxeN7P5MT6nAVcB17r7EOBa4P6GLmxmA4BDgVeiNt8IHAgcAfQGrm/oeHef6u757p7fr1+/OG83sTZu28XKTTvUfyIiHVq8fSgnu/v4qPWpZjbD3W8zs5tiHeDuJzR0MjN7GPhBsPokcF8j1z4HmO7uVVHnXhss7jKzB4kkuHarbkJIJRQR6cjiraHUmtk5ZpYRfM6J2uctuG4JkXnCACYASxopez71mruCWgsWaT86HZjfghjaTOGqLWRnGocMyg07FBGRpIm3hnIhMAX4U7D+PnCRmeUAV7fgulcAU8wsC6gEJgOYWT5wpbtPCtaHAUOAf9c7fpqZ9SMyhHkOcGULYmgzs1eWcfDAXLpkZ4YdiohI0sSVUNx9GbFnHAZ4p7kXdfd3gMNjbC8g8mbIuvUVRPpt6peb0NxrhmV3dS1zi8u46Kh9ww5FRCSp4h3lNTgY0bXBzNab2dNmNjjZwXUEC9duZVd1rZ4/EZEOL94+lAeBZ4GBRGoMzwXbpAmFekOjiKSJeBNKP3d/0N2rg89DQDhjcFNM4aotDOqZw4BcPdAoIh1bvAllo5ldZGaZweciYFMyA+soZq/cwmFDVTsRkY4v3oTyX0SeB1kHrCUyFYre4tiEkrKdlJRX6vkTEUkLcSUUd1/l7qe6ez937+/upwNnJjm2lFc3IaQSioikg5ZMX1+n0dcDC8xaWUaX7AwOGtCj6cIiIimuNQlFsxw2oXDVFsYM7kl2Zmv+M4uIpIbW/KVryZQraaOyqoYFJeV6/kRE0kajT8qbWQWxE4cBGgfbiHlryqmqcfWfiEjaaDShuHv3tgqko/n0gUYNGRaRNKHG/SSZtXILw/rsRZ9uDb77S0SkQ1FCSQJ3Z9aqLYxTc5eIpBEllCRYvXknG7ftVv+JiKQVJZQkKFy1GUAjvEQkrSihJEHhyi1065zFqL01pkFE0ocSShLMWlnGYUN7kpmhZz9FJH0ooSTYtl3VLFq3lcPU3CUiaUYJJcHmri6j1jUhpIikn9ASipnlmdn7ZjbPzJ4zs5gzKJrZyWa22MyWmtkNUduHm9lMM1tiZk+YWae2i75hs1ZuwQzGDtEDjSKSXsKsodwH3ODuhwLTgevqFzCzTOCPwNeA0cD5ZjY62P1r4G53HwlsAS5vk6ibULhqCyP7dyM3JzvsUERE2lSYCeUA4O1g+TXgrBhljgSWuvsyd98NPA6cZmYGTACeCsr9FTg9yfE2qbbWmb2qTM1dIpKWwkwo84FTg+VvAkNilBkErI5aLw629QHK3L263vZQLdu4jfKdVeqQF5G01OjkkK1lZq8D+8TY9RMirxX+vZn9DHgW2B3rFDG2eSPbY8UwGZgMMHTo0DiibrnPJoRUQhGR9JPUhOLuJzRR5CQAMxsFnBJjfzF71lwGAyXARqCnmWUFtZS67bFimApMBcjPz0/qO1zmrC6nR5cs9uvXNZmXERFpl8Ic5dU/+M4AbgbuiVHsQ2BkMKKrE3Ae8Ky7O/AmcHZQ7lLgn8mPunHz1pQxZnBPIl08IiLpJcw+lPPN7GNgEZHaxYMAZjbQzF4ECGofVwOvAAuBv7v7R8Hx1wM/NLOlRPpU7m/j+PdQWVXDorUVjBmcG2YYIiKhSWqTV2PcfQowJcb2EmBi1PqLwIsxyi0jMgqsXVi4divVtc6YwXr+RETSk56UT5Ci4nIA8oaohiIi6UkJJUHmFpfRt1tn9unRJexQRERCoYSSIEXF5eQNzlWHvIikLSWUBNi2q5pPSrep/0RE0poSSgLMKy7HHcao/0RE0pgSSgIUFZcBkKcaioikMSWUBCgqLmdwrxx6d20XM+iLiIRCCSUB5haXqXYiImlPCaWVNm/fTfGWnXpCXkTSnhJKK9X1n2iEl4ikOyWUVioqLscMDhkU8w3GIiJpQwmllYqKyxjRtyvdu+iVvyKS3pRQWsHdmVtcrg55ERGUUFpl3dZKSit2qUNeRAQllFaZuzoyw/CYIaqhiIgoobRCUXEZWRnG6AHqkBcRUUJphaLicg7YpztdsjPDDkVEJHRKKC3k7hQVl+n5ExGRgBJKC63YtIOtldXkqUNeRARQQmkxPSEvIrKnUBKKmeWZ2ftmNs/MnjOzz/Vqm9kQM3vTzBaa2Udm9oOofbeY2RozmxN8JrbtHUT6TzpnZTBq725tfWkRkXYprBrKfcAN7n4oMB24LkaZauBH7n4QcBTwXTMbHbX/bncfG3xeTH7IeyoqLuPggT3IylQlT0QEwksoBwBvB8uvAWfVL+Dua919VrBcASwEBrVZhI2orqll/pqtau4SEYkSVkKZD5waLH8TGNJYYTMbBhwGzIzafLWZFZnZA2bWKxlBNmRp6TZ2VtWQp1f+ioh8KmkJxcxeN7P5MT6nAf9FpAmrEOgO7G7kPN2Ap4Fr3H1rsPnPwH7AWGAt8JtGjp9sZgVmVlBaWpqQeyuqe0JeNRQRkU9lJevE7n5CE0VOAjCzUcApsQqYWTaRZDLN3Z+JOvf6qDL3As83EsdUYCpAfn6+xxt/Y+YWl9G9cxbD+3RNxOlERDqEsEZ59Q++M4CbgXtilDHgfmChu/+23r4BUatnEGlCazNFxeUcOjiXjAxry8uKiLRrYfWhnG9mHwOLgBLgQQAzG2hmdSO2vgBcDEyIMTz4jmDIcRHwZeDatgp8V3UNi9apQ15EpL6kNXk1xt2nAFNibC8BJgbL7wAxqwDufnFSA2zEwrUVVNW4npAXEalHD1E007y6J+Q1Zb2IyB6UUJppbnE5fbt1YmBul7BDERFpV5RQmqluhuHImAEREamjhNIM23dVs3TDNg4dpP4TEZH6lFCaYf6acmodPSEvIhKDEkozFBXrCXkRkYYooTTD3OIyBvXMoW+3zmGHIiLS7iihNENRcTlj9PyJiEhMSihx2rJ9N6s271Bzl4hIA5RQ4lS0JtJ/oifkRURiU0KJU9HqyBPyhyihiIjEpIQSp6I15Yzo15UeXbLDDkVEpF1SQolTUXEZeeo/ERFpkBJKHNZvrWT91l16Ql5EpBFKKHGYG/Sf6Al5EZGGKaHEoai4nMwMY/QAJRTVKKeQAAAIK0lEQVQRkYYoocRhSO8czh43mJxOmWGHIiLSboXyxsZUc+4RQzn3iKFhhyEi0q6phiIiIgmhhCIiIgkRSkIxszwze9/M5pnZc2bWo4FyK4Iyc8ysIGp7bzN7zcyWBN+92i56ERGJJawayn3ADe5+KDAduK6Rsl9297Hunh+17QbgX+4+EvhXsC4iIiEKK6EcALwdLL8GnNXM408D/hos/xU4PUFxiYhIC4WVUOYDpwbL3wSGNFDOgVfNrNDMJkdt39vd1wIE3/0bupCZTTazAjMrKC0tTUDoIiISS9KGDZvZ68A+MXb9BPgv4Pdm9jPgWWB3A6f5gruXmFl/4DUzW+TubzdQNiZ3nwpMBcjPz/fmHCsiIvFLWkJx9xOaKHISgJmNAk5p4BwlwfcGM5sOHEmkqWy9mQ1w97VmNgDYkLjIRUSkJUJ5sNHM+gdJIgO4GbgnRpmuQIa7VwTLJwG3BbufBS4Fbg++/xnPdQsLCzea2coWht0X2NjCY9sb3Uv701HuA3Qv7VVr7mXfeAqZe9u3ApnZD4DvBqvPADe6u5vZQOA+d59oZiOIjACDSOJ71N1/ERzfB/g7MBRYBXzT3TcnOeaCeiPNUpbupf3pKPcBupf2qi3uJZQairtPAabE2F4CTAyWlwF5DRy/CfhKMmMUEZHm0ZPyIiKSEEoo8ZsadgAJpHtpfzrKfYDupb1K+r2E0ociIiIdj2ooIiKSEEoocTCzk81ssZktNbOUnjesoQk32zsze8DMNpjZ/KhtKTlJaAP3couZrQl+lzlmNjHMGONlZkPM7E0zW2hmHwUjOFPut2nkPlLudzGzLmb2gZnNDe7l1mD7cDObGfwmT5hZp4RfW01ejTOzTOBj4ESgGPgQON/dF4QaWAuZ2Qog391Tamy9mX0R2AY87O6HBNvuADa7++1Bou/l7teHGWc8GriXW4Bt7n5XmLE1V/Bg8QB3n2Vm3YFCInPrXUYK/TaN3Mc5pNjvYmYGdHX3bWaWDbwD/AD4IfCMuz9uZvcAc939z4m8tmooTTsSWOruy9x9N/A4kckppQ0FU+7Uf9YoJScJbeBeUpK7r3X3WcFyBbAQGESK/TaN3EfK8YhtwWp28HFgAvBUsD0pv4kSStMGAauj1otJ0f+hBRqacDMVxT1JaIq42syKgiaxdt1EFIuZDQMOA2aSwr9NvfuAFPxdzCzTzOYQmZbqNeAToMzdq4MiSfk7poTSNIuxLZXbCb/g7uOArwHfDZpfJHx/BvYDxgJrgd+EG07zmFk34GngGnffGnY8LRXjPlLyd3H3GncfCwwm0spyUKxiib6uEkrTitlzev3BQElIsbRa9ISbRKa2OTLciFplfdD2XdcGnrKThLr7+uCPQC1wLyn0uwTt9E8D09z9mWBzyv02se4jlX8XAHcvA94CjgJ6mlnd7ChJ+TumhNK0D4GRwQiJTsB5RCanTDlm1jXocKybfPMkIu+mSVV1k4RCMyYJbY/q/vgGziBFfpegA/h+YKG7/zZqV0r9Ng3dRyr+LmbWz8x6Bss5wAlE+oTeBM4OiiXlN9EorzgEQwV/B2QCD9RNUplqGptws70zs8eA44nMmLoe+DnwD9p4ktBEaOBejifSrOLACuDbdX0Q7ZmZHQv8B5gH1AabbyLS/5Ayv00j93E+Kfa7mNkYIp3umUQqDX9399uC//8/DvQGZgMXufuuhF5bCUVERBJBTV4iIpIQSigiIpIQSigiIpIQSigiIpIQSigiIpIQSiiStsxsW/A9zMwuSPC5b6q3/l4iz1/v3NeY2V7JOr9IvJRQRGAY0KyEEsxC3Zg9Eoq7H9PMmJrjGkAJRUKnhCICtwPHBe+7uDaYWO9OM/swmBTw2wBmdnzwzoxHiTwAh5n9I5ho86O6yTbN7HYgJzjftGBbXW3oieh3apjZQ2Z2VkPXjBbMdPBC8J6L+WZ2rpl9HxgIvGlmbwblTjKz981slpk9GcxPVfcunF8H78r4wMz2D7Z/MzjfXDN7O2n/laXjc3d99EnLD5H3XEDkKfXno7ZPBm4OljsDBcDwoNx2YHhU2d7Bdw6RaTn6RJ87xrXOAP4aLHciMpN1TkPXrHeOs4B7o9Zzg+8VQN9guS/wNpH3YQBcD/wsqtxPguVL6u6ZSHIcFCz3DPt30Sd1P6qhiHzeScAlwfTfM4E+wMhg3wfuvjyq7PfNbC4wg8gkoiNp3EvABDPrTGTG57fdfWcT16wzDzghqGUc5+7lMc5/FDAaeDc416XAvlH7H4v6PjpYfhd4yMyuIDJdh0iLZDVdRCTtGPA9d39lj41mxxOpoUSvnwAc7e47zOwtoEtjJ3b3yqDcV4Fz+ewPfMxr1jv2YzM7HJgI/MrMXnX322LE/pq7n9/Qaeovu/uVZjYeOAWYY2Zj3X1TY/chEotqKCJQAXSPWn8FuCqYzhwzGxXMzlxfLrAlSCYHEqkd1KmqOz6Gx4FvAccF14rrmmY2ENjh7o8AdwHjYsQ/A/hCVP/IXmY2Kuo050Z9vx+U2c/dZ7r7z4CN7Pm6BpG4qYYiAkVAddB09RAwhcjIr1nBtOalxH5d6svAlWZWBCwm8se8zlSgyMxmufuF9Y57FXgYeNYjr5UGuC+Oax4K3GlmtUAVcFXUtV4ys7Xu/mUzuwx4LGhWA7gZ+DhY7mxmM4n8Y7KuFnOnmY0kUrv5FzA35n8lkSZotmGRNGFmK4B8d98YdizSManJS0REEkI1FBERSQjVUEREJCGUUEREJCGUUEREJCGUUEREJCGUUEREJCGUUEREJCH+H7OWtFexyyQjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gmm_salinas = GMM(potential_features=pCs_salinas, n_coms=6, em_steps=50, n_features=2)\n",
    "gmm_salinas.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(Y, gmm):\n",
    "    Y = Y[:, :gmm.n_features]\n",
    "    labels = np.zeros((Y.shape[0], gmm.n_coms))\n",
    "    for c in range(gmm.n_coms):\n",
    "        labels [:,c] = gmm.pi[c] * multivariate_normal(gmm.mu[c,:], gmm.cov[c]).pdf(Y)\n",
    "    labels  = labels.argmax(1)\n",
    "    return labels\n",
    "    \n",
    "\n",
    "def compute_confusion_matrix(true, pred):\n",
    "    K = len(np.unique(true))\n",
    "    print(np.unique(true))\n",
    "    result = np.zeros((K, K))\n",
    "\n",
    "    for i in range(len(true)):\n",
    "        result[true[i]][pred[i]] += 1\n",
    "\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEXCAYAAACJYMEPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH05JREFUeJzt3X+sXOV95/HPJ/4BGALG2CDHJgtRvAkhy4/sFSXJqmIhcQJFAXXDQkJbp2Llf7otTtMttH9skmorJWpV0pVWqFZI60r5AaUkIERCWBf2h9S4MYE6ASfBpQQcHGwIJvzogg3f/WPOtc81d7hn5nnOmXPOvF8S8sx47p3neqyHx897zjmOCAEAAGA8b5r0AAAAALqMxRQAAEACFlMAAAAJWEwBAAAkYDEFAACQgMUUAABAAhZTAAAACZIWU7Y/bPtHtnfZvj7XoACgCcxhAHLwuCfttL1I0o8lfVDSbknflfSxiHg43/AAoB7MYQByWZzwtedJ2hURj0qS7a9JukzS0Ilo6ZuOiWMWvznhJTGuOHBw0kM4ZNkZkx5Bc17amef7nHTmK3m+UcOe/unLev7ZA570OIYYaQ5buWJRnHbqkgaHh1k/3rFs0kM45OVTjz10+5Tjn5vgSOq3/6GUJcJh5T+zrnnlid1PR8SqhZ6X8ie1RtITpfu7Jf3SkU+yvVHSRkk6etFxet+qKxNeEuM6uOdnkx7CIWd9pa3/b81vx3vyXK7p6tt2Z/k+Tfvsr35/0kN4IwvOYeX5661rFusf7j61udHhkA+95ZxJD+GQXZ86/9DtTeu/OcGR1O/OM0/M8n3Kf2Zd89im3/tJleelfGZqvv8jvu7/HBGxOSJmImJm6ZuOSXg5AMhqwTmsPH+tOmlRQ8MC0DUpO1O7JZX/mbZW0pNpw0GqVu1AfY8dqBRX/7Cbu1EdwhzWMq3agbph2A5U/3ajcu1AlZX//KZBys7UdyWts3267aWSrpJ0R55hAUDtmMMAZDH2zlREHLT9nyXdLWmRpC9FxEPZRgYANWIOA5BL0kf1I+IuSXdlGgvGNOm0R85LQ86bHOawyZt02pumHEXOqw9nQAcAAEjAYgoAACBBnjNyoXGkvXqR84D6tCnt9fFcUeS85rEzBQAAkIDFFAAAQAIyX4dMOu31XbZLv5DzgNeZdNrru2yXfiHnjYWdKQAAgAQspgAAABKQ+VBZH4/gI+0B06GP+Yq01x7sTAEAACRgMQUAAJCAzNdyHMEHoKvadARfH0/OifZgZwoAACABiykAAIAEZL4WalPa4wg+AKNoU9q79KFnJz2E7Oq47h7SsTMFAACQgMUUAABAAjJfS7Qp7fVRHWmPE3UCA21Ke31UR9rjRJ15sTMFAACQgMUUAABAAjJfw9qa8/p01B5H6wH1aGvO69NRexyt100L7kzZ/pLtvbZ/UHpshe17bD9S/Mq7D6CVmMMA1K1K5vsrSR8+4rHrJW2NiHWSthb3AaCN/krMYQBqtGDmi4j/bfu0Ix6+TNIFxe0tku6TdF3GcXUeOa9+TeU8jtrrNuaw0ZHz6tdUzuOovWaM+wH0UyJijyQVv5487Im2N9rebnv7K6/9y5gvBwBZVZrDyvPXvmdebXSAALqj9g+gR8RmSZsl6YSlJ/fuk8Ft3YGS+rMLxQ4UJqU8f82cfXTv5q+27kBJ/dmFYgdqOoy7M/WU7dWSVPy6N9+QAKB2zGEAshl3MXWHpA3F7Q2Sbs8zHABoBHMYgGwWzHy2v6rBBzVX2t4t6dOSPifpFtvXSHpc0hV1DrIN2pzzykh7CyPn1e/L71wrSXomfjThkfR/DvvQr/7GoduX3vS/Dt3+7RN/cvg5Lc55ZaS9hZHz6vf2T37n0O3HKn5NlaP5Pjbkty6q+BoAMDHMYQDqxuVkAAAAEnA5mTfQ5rRHzquGpFev2ZyHybn7tr+e9/E2pz1yXjUkvXqVc14qdqYAAAASsJgCAABIQOYTOa8J5LxuI+e1V5tzXjlTbVr/zQmOJA05r9ty5rxh2JkCAABIwGIKAAAgwVRlPnJes0h73Ubaa5c257zhR+eR9oYh7dWribRXxs4UAABAAhZTAAAACXqb+dqc9Gb1Je2R87qNnNc+bU56szjxZjXkvHo1nfOGYWcKAAAgAYspAACABJ3PfOS85pDzuq3JnPcnjw223j9+6QuNvWYXkfOaQ87rtiZz3tMb33v4zl/cWulr2JkCAABIwGIKAAAgQSczX1vTHjmvGnJe/ZpKerM5D9W1Ne1xHb1qyHn1ayrpzcl5idiZAgAASMBiCgAAIEGrMx85rzl1pD1yXrPqTHvkvNF1L+eR9srIec2qM+3lzHnDLLgzZftU2/fa3mn7IdvXFo+vsH2P7UeKX+sN1QAwIuYvAE2okvkOSvpURJwh6XxJv2X7XZKul7Q1ItZJ2lrcB4A2Yf4CULsFM19E7JG0p7j9vO2dktZIukzSBcXTtki6T9J14wyCnNccjtTrPnJedU3MX23NecNPtknOG4a0V7+u57xhRvoAuu3TJJ0raZukU4qJanbCOnnI12y0vd329lde+5e00QLAmFLnr33PvNrUUAF0TOXFlO3jJP2tpE0R8YuqXxcRmyNiJiJmlr7pmHHGCABJcsxfq05aVN8AAXRapaP5bC/RYCL6ckTcVjz8lO3VEbHH9mpJexf6PnHgYGuT3qw+pb06kx45r36ceDOPXPPXj3csa23SmzU3U3U350n1Jr225bxhJ0n9wrcvbngk+XTxxJspqhzNZ0k3SdoZEX9W+q07JG0obm+QdHv+4QHA+Ji/ADShys7U+yX9uqTv236weOwPJX1O0i22r5H0uKQr6hkiAIyN+QtA7aoczfd/JQ1rXxflHU5z+pTzynKnPXJe/ch59enr/NXHI/Wk/GmvbTlvGHLewtqS84bhcjIAAAAJWEwBAAAkaPW1+XLpY9Ij53Vb3WlvGpNeXw1Pet01rTmvrHwEX9cyX91pr+1Jbz7sTAEAACRgMQUAAJCgV5mPnLcwcl6zyHmoaniq6u4RetOU84adeLOsnPO6kPbIedWxMwUAAJCAxRQAAECCzme+vqS9Oq6jR9JrDjlvrv9y2iDH7I6tEx5Ju5WzVZVM1FZ1XEevzUlvPl3IdsOQ86Qll+87dPuES3aN/PXsTAEAACRgMQUAAJCgM5mvLzmvjLTXbaS9wzkPb6yP19KbprQ3LMGS9obrQtpbufnvD9/ZnPa92JkCAABIwGIKAAAgQesyHzmvGnJec8h55LyqOPFmNeS85pDzjsh5NWFnCgAAIAGLKQAAgAQTy3zkvGrIec0h581F2huOo/OqaWvO6yNy3lxNpL0ydqYAAAASsJgCAABI0GjmW3aGdNZX+pX3cqU9cl5zyHlz5cp5dz/54OseO+9DL2X53m2w/MyDuvRvhuW9bsqV9rqQ86pc+7ALR+2R8+bKlfOeu+vt8/9Gxb8SC+5M2T7a9j/Y/kfbD9n+bPH46ba32X7E9s22l1YfNgDUj/kLQBOqZL6XJV0YEWdLOkfSh22fL+nzkm6IiHWSnpV0TX3DBICxMH8BqN2CmS8iQtILxd0lxX8h6UJJHy8e3yLpM5JuzD/EdiDndRtpr96c11bMXwN9z3l9SXjDkPYayHmJKn0A3fYi2w9K2ivpHkn/JGl/RBwsnrJb0pohX7vR9nbb21/a/3KOMQNAZbnmrxd+fqCZAQPonEqLqYh4NSLOkbRW0nmSzpjvaUO+dnNEzETEzLLlR40/UgAYQ67567gVS+ocJoAOG+lovojYb/s+SedLWm57cfGvu7WSnqxhfI0j53VfnUmvCzlPyn/CzS6lvWGmYf7qe86bBnUmvS7kPCn/CTfrSntlVY7mW2V7eXH7GEkfkLRT0r2SPlo8bYOk2+saJACMg/kLQBOq7EytlrTF9iINFl+3RMSdth+W9DXb/03SA5JuqnGc2XHpl+5jB4odqAp6OX9N66Vf+vRBc3agurkDNUyVo/l2SDp3nscf1eDzBwDQSsxfAJrA5WQAAAASNHo5mUnjw+XdRM6bK0fa62HO671p+3D5QkmPnNednFeWI+1NMucNw84UAABAAhZTAAAACXqb+XIkPXLe5JD2pvPyLxjIkfS6kvP6iLTX/su/5MbOFAAAQAIWUwAAAAk6n/nIed1WZ86TupH0yHnTa1pzXvlIvWFH5XXhaL06c57UjaQ3bTlvGHamAAAAErCYAgAASNDJzDdu2iPnTQ45L/919CTSXheNm/a6kvOGnWyznO26kPDKyHn5r6MndT/tlbEzBQAAkIDFFAAAQIJWZ76UI/VIepNBziPnYSDlSL2uJL2+IeeR8yTphD9aNvLXsDMFAACQgMUUAABAglZkvhwn3pRIe5NSR9rrQs4r48Sbc33oLedIkn4cz0x4JPXLceJNqRtpjyP1qulCziub1hNvDs1539kx8vdiZwoAACABiykAAIAEjWa+l3ZyLb0u40g9ct6RZnPeNNj/0OKpuZbesGvndS3hlXGk3vTmvKHGyHnDVN6Zsr3I9gO27yzun257m+1HbN9se2m2UQFARsxfAOo0Sua7VtLO0v3PS7ohItZJelbSNTkHBgAZMX8BqE2lzGd7raRfkfTHkn7XtiVdKOnjxVO2SPqMpBtzDo6cNxnkPHLekbqc8yY1f3Uh5/XRtOa8JZfvO3T7hEt2ZfmeXct5OY/OG1XVnakvSPp9Sa8V90+StD8iDhb3d0taM98X2t5oe7vt7Qf0ctJgAWAMzF8AarXgYsr2pZL2RsT95Yfneeq8nyyPiM0RMRMRM0t01JjDBIDRMX8BaEKVzPd+SR+xfYmkoyUdr8G/9JbbXlz8626tpCfHHQQ5bzLIeVxH7410Oe2V1D5/9THndeGovWnNeWVzjs7bnOd7di3tDdVA2itbcGcqIv4gItZGxGmSrpL0dxFxtaR7JX20eNoGSbfXNkoAGAPzF4AmpJy08zoNPsy5S4PPINyUZ0gAUDvmLwDZjHTSzoi4T9J9xe1HJZ03ytefdOYruvo2kt4kTet19Mh5w+XKeUvuW/26x7xxSZbvnUPq/PXyqcdq16fakfSqXCOvj6b1Onq5TrZZ1uWcN+eovVw57/yz5n/872+t9OVcTgYAACABiykAAIAEjV6bD82q82i9LqQ9KU/eI+fNNV/OA3Kr82i9LqQ9KU/em1TOO/CNVYdur9zx4rzPee6/vjTv4/OefLPunJeInSkAAIAELKYAAAASkPl6gJxHzjtS7hNukvYmb9gRfF1Hzut2zhumfK3A5y6f/zmNXUuvprRXxs4UAABAAhZTAAAACch8HULOI+cdiZzXfwvlva6cqJOc18+cl6SDOW8YdqYAAAASsJgCAABIQOZrudxpb5pyXhlpby5yXnd0JePNJ3fam6acV9bFtFfrkXoTzHnDsDMFAACQgMUUAABAAjJfS9RxpF4Xkl7unCf1J+lxLT10RR1H6nUh6eXOeVI3k96sOWmv5dfSy42dKQAAgAQspgAAABKQ+RpW54k3pelKe+S8uch502nYST3rOAqwzhNvStOV9rqW82q/jl6Lct7TZx17+E7Ft5udKQAAgAQspgAAABKQ+WpCziPnHSn3dfQk0h7IeXWZ1pw3VO7r6EmtSnu7rjycMTetv+vw7b+o9vWVFlO2H5P0vKRXJR2MiBnbKyTdLOk0SY9J+o8R8Wy1lwWAZjB/AajbKJnv30fEORExU9y/XtLWiFgnaWtxHwDaiPkLQG1SMt9lki4obm+RdJ+k6xLH02mceDOfLqc9cl4nMH8dgRNv5tPltHfCJbvyf9MJ5Lw5R+RVEkmvV3VnKiR92/b9tjcWj50SEXskqfj15Pm+0PZG29ttb3/+2QNJgwWAMWSZv1594cWGhguga6ruTL0/Ip60fbKke2z/sOoLRMRmSZsl6fR3H5e29AOA0WWZv45666nMXwDmVWkxFRFPFr/utf11SedJesr26ojYY3u1pL01jrNVyHn5dDnnlXHizfZi/pqLnJdP13LeNJ14c+WOwzvJ5eT3iWvvmu/pyRbMfLaPtf3m2duS1kv6gaQ7JG0onrZB0u21jBAAxsT8BaAJVXamTpH0dduzz/9KRHzL9ncl3WL7GkmPS7qivmFOxrSeK4odqOHYgeqcqZ2/pvVcUexAzTVnN6rDO1Cjf6D8sP1n1F/oF1xMRcSjks6e5/FnJF1Ux6AAIAfmLwBN4HIyAAAACbicjKY355Vx6Ze5yHnoimnNeWVc+mWubOeKmtAHyhdKesOy3ab136xjOJWwMwUAAJCAxRQAAECCqc18nCsqDTlvONIe6sa5otKQ895Ai84V1SXsTAEAACRgMQUAAJBgKjJf7qTXhZwn5T/5ZpfTHjkPXZU76XUh50n5T77Z5bQ3zTlv9si9SR6pVwU7UwAAAAlYTAEAACToVeYj5+XR5Zx3pBx5j5yHJpDz8uhyzpNquJZeS0+8Kb3RyTfvyj2cQ77w7YuHvGZaRmRnCgAAIAGLKQAAgASdzHx1XkuvzWmPnDccOQ9dUee19Nqc9sh5c3U555UT3sodLy74/ElfS6+c9t5+80uHbu+6ctl8Tx8LO1MAAAAJWEwBAAAk6Ezm40i9NH1JerlPvknaQxM4Ui9Nl5Ne9pxXVnPaq3JE3rDnlNPepE+4Oef119fzGuxMAQAAJGAxBQAAkKB1ma+OI/W6kPTIeXOR89BFdRyp14WkN4mcd+Abq+Z9fMnl+7KMZVRdznmjmvTReW1UaWfK9nLbt9r+oe2dtt9re4Xte2w/Uvx6Yt2DBYBRMX8BqFvVzPfnkr4VEe+UdLaknZKul7Q1ItZJ2lrcB4C2Yf4CUKsFM5/t4yX9sqRPSFJEvCLpFduXSbqgeNoWSfdJuq7qC5Pz8uly2sud8ySSHg6ra/4i5+XTtSP1as150kSS3qjX0ZvmnDdMlZ2pt0naJ+kvbT9g+4u2j5V0SkTskaTi15Pn+2LbG21vt739+WcPZBs4AFSQbf569YWFz/QMYDpVWUwtlvQeSTdGxLmSXtQIW+IRsTkiZiJi5s0nLhlzmAAwlmzz16LjFv7XO4DpVOVovt2SdkfEtuL+rRpMRk/ZXh0Re2yvlrR3oW/0zENLp+bkm+S84XKlPXIeKsg2fx31xItTc/LNLuS8SR21ly3ttewIPaRZcGcqIn4m6Qnb7ygeukjSw5LukLSheGyDpNtrGSEAjIn5C0ATqp5n6rclfdn2UkmPSvpNDRZit9i+RtLjkq6oZ4gAkIT5C0CtKi2mIuJBSTPz/NZFeYczHDmve8h5aIM2zF/kvO454ZJdeb7RhHNelSP1yjhqbzxcTgYAACABiykAAIAErbs2X1tzXhnX0ZuLE28CA23NeWWTuI5em2XLeWUdSHvkvLzYmQIAAEjAYgoAACDBxDIfOa97yHnAADmve/qS88oJb+WOF+d9vJzwysh59WFnCgAAIAGLKQAAgASNZr61/+YF/cmd7ct7nHhzLnIe8HoHVx6rp/9D+/IeJ96cq+s5b74j8co5r3x715XLDt3etP6uegeGN8TOFAAAQAIWUwAAAAlad9LOOpHzhuM6ekC7kfOG68t19NBd7EwBAAAkYDEFAACQoLeZL3fSI+cNR9oD8sqd9Mh5b6AFaW/YtfRmT765/4zyUXuceLON2JkCAABIwGIKAAAgQa8yX46015ecJ+VPeuQ8VHXggj2SpIgDEx5Jd+RIe7ly3oFvrDp0e8nl+7J8z1FlT3oTPvHmkeZeP+/wbTLe5N155okjfw07UwAAAAlYTAEAACToZOYj581FzsOkzOY8VNemnNcGXc556LZxct4wC+5M2X6H7QdL//3C9ibbK2zfY/uR4td8owKADJi/ADRhwcVURPwoIs6JiHMk/VtJL0n6uqTrJW2NiHWSthb3AaA1mL8ANGHUzHeRpH+KiJ/YvkzSBcXjWyTdJ+m6fEPjxJtHIudhUnqS8xqdv7p84s3yEXy5juzrY84rH7W3cseL8z6nfNQeR+pNRs6cN8yoH0C/StJXi9unRMQeSSp+PXm+L7C90fZ229v3//y18UcKAGmS5q+D/2/+/1kCQOXFlO2lkj4i6W9GeYGI2BwRMxExs3wFBw8CaF6O+Wvx0QufOwjAdBol810s6XsR8VRx/ynbqyNij+3VkvaOO4jcOU/qftKblSvtkfQwqp6kvVm1zV+5c57Un6P1sqW9FiQ9dEsTaa9slK2ij+nwFrkk3SFpQ3F7g6Tbcw0KADJj/gJQm0qLKdvLJH1Q0m2lhz8n6YO2Hyl+73P5hwcAaZi/ANTNEbHwszI53ivil3xR1u9JzhuOtIcq6sx522KrfhE/d20v0KA65i9y3huYcNob/fp6h3HUXnPqznn/M269PyJmFnoenwgHAABI0JnLybADNRw7UKiiZx8o75Qu70Cd8EfLDt/5zo78L1DjDlSV80ANez7nh2qXpj9QPip2pgAAABKwmAIAAEjQusxHzhuOnIdxkPea0+WcV9bHD5Sjm9qe92axMwUAAJCAxRQAAECCVmS+Lqc9ch7agpw3GV1Oe9OU84YdqVe2af1dTQ0HR+hKzhuGnSkAAIAELKYAAAASNJr5/vVZL+nuu7ub9GblSnvkPIyDnDcZr647Ss/99+4mvVnZ0l5Lcx7ares5bxh2pgAAABKwmAIAAEjQiqP52ooj9TBJ5Dyk6PuRehyd1259zXnDsDMFAACQgMUUAABAAjKfyHloD9IeRlVHznt643vnfXzljhezvxb6Y9rSXhk7UwAAAAlYTAEAACSYqsxHzkNbkPMwqiaPzptUzisfoVc2e7TepvXfbHI4GGKac94wlXambH/S9kO2f2D7q7aPtn267W22H7F9s+2ldQ8WAEbF/AWgbgsupmyvkfQ7kmYi4t2SFkm6StLnJd0QEeskPSvpmjoHCgCjYv4C0ISqmW+xpGNsH5C0TNIeSRdK+njx+1skfUbSjbkHOK7cSY+ch3GQ81qhc/NX9qTXopNtojvIedUtuDMVET+V9KeSHtdgEnpO0v2S9kfEweJpuyWtme/rbW+0vd329n3PvJpn1ABQQc756+BzLzUxZAAdVCXznSjpMkmnS3qLpGMlXTzPU+c9n39EbI6ImYiYWXXSopSxAsBIcs5fi09YVt9AAXRalcz3AUn/HBH7JMn2bZLeJ2m57cXFv+7WSnqyvmEOV8cRerNIexgHaa9VWj1/1XKE3qyWpb1hR+p94lqun9cmpL3xVDma73FJ59teZtuSLpL0sKR7JX20eM4GSbfXM0QAGBvzF4DaVfnM1DZJt0r6nqTvF1+zWdJ1kn7X9i5JJ0m6qcZxAsDImL8ANKHS0XwR8WlJnz7i4UclnZd9RBVwpB7agJzXDW2bvzhSD21AzsuLy8kAAAAkYDEFAACQwBHzHhFcz4vZ+yS9KOnpxl50slaKn7WP+Fmr+1cRsSrXYCapmL9+It7/vuJn7adG5rBGF1OSZHt7RMw0+qITws/aT/ys022a/kz4WfuJnzU/Mh8AAEACFlMAAAAJJrGY2jyB15wUftZ+4medbtP0Z8LP2k/8rJk1/pkpAACAPiHzAQAAJGh0MWX7w7Z/ZHuX7eubfO262T7V9r22d9p+yPa1xeMrbN9j+5Hi116cdtb2ItsP2L6zuH+67W3Fz3mz7aWTHmMutpfbvtX2D4v39709fl8/Wfz9/YHtr9o+us/v7SiYv/rz91yanjmM+auZ97WxxZTtRZL+h6SLJb1L0sdsv6up12/AQUmfiogzJJ0v6beKn+96SVsjYp2krcX9PrhW0s7S/c9LuqH4OZ+VdM1ERlWPP5f0rYh4p6SzNfi5e/e+2l4j6XckzUTEuyUtknSV+v3eVsL81Z+/5yXTMocxfzXwvja5M3WepF0R8WhEvCLpa5Iua/D1axUReyLie8Xt5zX4C7tGg59xS/G0LZIun8wI87G9VtKvSPpicd+SLtTggrJST35OSbJ9vKRfVnEh3Ih4JSL2q4fva2GxpGNsL5a0TNIe9fS9HRHzV4/e+2mZw5i/mpu/mlxMrZH0ROn+7uKx3rF9mqRzJW2TdEpE7JEGE5akkyc3smy+IOn3Jb1W3D9J0v6IOFjc79N7+zZJ+yT9ZZEEvmj7WPXwfY2In0r6U0mPazAJPSfpfvX3vR0F81dP/p4XpmUOY/5q6H1tcjHleR7r3aGEto+T9LeSNkXELyY9ntxsXyppb0TcX354nqf25b1dLOk9km6MiHM1uBxS57fE51N8buIySadLeoukYzXIWkfqy3s7ij7/HT+k7/OXNHVzGPPX69Xyvja5mNot6dTS/bWSnmzw9Wtne4kGE9GXI+K24uGnbK8ufn+1pL2TGl8m75f0EduPaZA6LtTgX3nLi61VqV/v7W5JuyNiW3H/Vg0mp769r5L0AUn/HBH7IuKApNskvU/9fW9HwfzVn7/n0zSHMX819L42uZj6rqR1xSfrl2rwwbA7Gnz9WhXN/SZJOyPiz0q/dYekDcXtDZJub3psOUXEH0TE2og4TYP38O8i4mpJ90r6aPG0zv+csyLiZ5KesP2O4qGLJD2snr2vhcclnW97WfH3efZn7eV7OyLmr56899M0hzF/NTd/NXrSTtuXaPAvgEWSvhQRf9zYi9fM9r+T9H8kfV+HO/wfavC5g1skvVWDN/uKiPj5RAaZme0LJP1eRFxq+20a/CtvhaQHJP1aRLw8yfHlYvscDT6oulTSo5J+U4N/iPTufbX9WUlXanB01wOS/pMGnzHo5Xs7Cuav/vw9nzUNcxjzVzPzF2dABwAASMAZ0AEAABKwmAIAAEjAYgoAACABiykAAIAELKYAAAASsJgCAABIwGIKAAAgAYspAACABP8fiWDgYH2ic7wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels = predict(pCs_salinas, gmm_salinas)\n",
    "\n",
    "count, loc, size = 0, 0, salinas_gt_flat.shape[0]\n",
    "final = np.zeros(size)\n",
    "while count<labels.shape[0]:\n",
    "    if (salinas_gt_flat[loc]!=0):\n",
    "        final[loc]=labels[count]+1\n",
    "        count+=1\n",
    "    loc+=1\n",
    "\n",
    "final = np.reshape(final, (n, m))\n",
    "\n",
    "f, ax = plt.subplots(1, 2, figsize=(10, 10))\n",
    "ax[0].imshow(salinas_gt)\n",
    "ax[1].imshow(final)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/joshua/anaconda3/lib/python3.7/site-packages/sklearn/metrics/cluster/supervised.py:732: FutureWarning: The behavior of AMI will change in version 0.22. To match the behavior of 'v_measure_score', AMI will use average_method='arithmetic' by default.\n",
      "  FutureWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8530955317687071\n"
     ]
    }
   ],
   "source": [
    "def randIndex(gt_labels, pred_labels):\n",
    "    return adjusted_mutual_info_score(np.reshape(gt_labels, -1), np.reshape(pred_labels, -1), average_method='warn')\n",
    "\n",
    "print(randIndex(salinas_gt, final))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}