MLMI_Exercise02_KMeans_PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Group A\n",
    "This code is for understanding the **k-means** algorithm. \n",
    "Follow the instructions mentioned before every function and run the codes for \n",
    "* 1, 5, 15 iterations \n",
    "* For 5 different intialization of centroids \n",
    "* Number of clusters : 2, 3, 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADaCAMAAABqzqVhAAAB41BMVEX/////zJn//wDM////mcwA/wD/AAAAAAD/nND/0Jz/0p3G+fl6Wzt/enqZUXXs7+7CwsI+Pj719fXT///k5ADa2tpWVlYlJSW5kWhxWkT5+vxxRFqTk5NacXHLy8tZZmaIk5OxZYuYm5nk4+OLi4tjYz1kZEhKShWlpZJ6iooA+AAA8AARERFLS0vmirj/oddeXl6+AAAA5QDPzwCCoqIA1wAAVACkpKSbwsJycgAAtgDu7gDa2gBEVVXWgKueAAAAjgDtAAAAzAC86+sAagDIeKDKAACiYoKvcpCvkHKigmJgYAC1bZGVWXcAUQCr1tZqamo1NTV7ewAAvwBsAACwsrHHqKiLYmIAqgAAYACKinRaNkhSUgBWVioAfACrAABHKzmQc1ZmPVIAQgAAhgDPpXyZmQC1tQBlf398SmNvi4sfHwA9TU0AcwB4eGB3AABHAAAALwCDAABVVQCQkAC9vQBWRDSti4smFx4AmwAnSyc9nD09Zz1rYWtmFBSKLy8AFwBkQ0NLPEtvsG8udy52nHYyAB6Fk4XFzsXOxM5XuVcnMicfBwA6OgC+vnN6enJ9fT8/MiWtaGg5AACSGhprMTEsLAChoU96ejFXalc4cziYTEyVypWksKQIIwhMGhooYygcHpgzAAAU4klEQVR4nOVdi18iV5ZWHlI26/ISVBDoiEk6oKW8HwoqZUlGCvFFAsEo4BMi07atTrtJemfz2Elv0tuT7ZnMpmc2u3/qnlvFU5CHQBXo9/vdoupwkftxzrn3u9d69PU1hMeffPrpJ6j84dM/0OVTpaGxT/YUnHtLN2H+q4TrVrUdZ0sDVbAwyXW72gwn0DQ/fboysFTkuwLl7/fLoxEz0Ho7sPR2ZQH2GKoLC7BZest109qKT4DayhnQWjo7W1rY3l4yn2+f0aH8lOumtROSc9qBZ2fmgQXz0gK4cu984Bz5c8Ac4bpxbQTDc2BpZdu8Zx7Ye/pyAV4YnitSrhvXRtA8zQuoLJjN5wN7JTzv0xiq/i/kzZfne2dL5u297b2X29Admen8XLhP/uyLrNBhu7KEOlt4oUcXuts947ppbYXyZTWVANgb4rpp7YVhryrNpXtGEzz6lyrC7755E8HwdnuBwdPc6/nLe0gTIRIxGCJSicypM+giOsO96mkrYZjiugVsAObW6zqD/D6pg6qIjMyKRmRD957ntAhBxnUzOg71LOJ5z7sghEnEU851KzoPxPOejptlcALP+zS3vhVO0ca9720RZkUqrpvACgwPobftQzzHuW5COyHVSQw6VJQSpc4g0dEFmeRO2nTTXKUm1xQawfjOh++//+Hn739eKHCYLzVMJeZfeuCfEtILXstwc02iAfww3zLNeRHXJOpDt9m6O/+7B0bZr4dbpjm3zjWJ+pD/Ld/a4S+GeXM7payH3ZWc5n85qrC5nVyzqI9vrAVWok3gOVf0EuKJDoeZI2bPelxBc/gJ1yQaQLDY3GDwmLczH+RZNy9fu3/ZdLvnvnAHj+d23Dvzm68353hzH27uXD7/sCLOL3pANDmtRZ471p25nXk38LRu8jYvecH55+DgoyPe5ebFJeKD9uYvK/wp6/4pauSXkvDb4R27n9M8L494m8PAE+IWvOl2H10gJ7qPedbn1gqeR1yTaACRyzKevJ0vIG4v3Xme4M8g+HP4eHM+58+Livyc+4JrEvWhLul26N51eGd4c8d9BGTAg5vDO+6gdc7tDs4d0UkZdAeHK+L2sge0rbNijADepfJoDv0ORUM15dQDna36gyrtbhKXPTBzq3Rn89jhmkQDeFJV8s1VM+bsFe9d9sDYqa4+IXMXe5pLa4kdFMVmxSf+3AMnib2p4s75YYYnqB/rHA8EEs8KtebmaZ4XFXp3uAeWd5VV3Lm56T52W4No7Nw5Cl7uBOeCF0Hr8S8XELDBo0pZH+yBsHVWjhIgEmAYtbrBqcOvwZ9uq/WCZw0e0/2VyP38Znpaf+CaRH2o/7PSnVak6nM8rRevre7L4+DFxfExnbE7vAqHbiq5ZlEfb6yVPOee84YhXN1W3mvrEe/4yH0JEnf+6JjWevAb7JSLobkhNdcs6uPPlTR5vONgcP5o3grSb24zCD4FxRcMDjPdLqTzfHlKX/TAKt+rKu7k5YReYVJdtFSt+7r7BxXDzm2tbxxHPbD6Jalc/mgaU93fC0met07zuAfGzsk2uPNzrknUh/rN/HBDqFHtqBcWbWeffPbBZ09+T5cnTz5A5fe06bNykyhnrlLzDdckGoSaKWp1blddZupjDmflt9e8R5BPd/+SZTtgmHoYPCU9sATdDqhnHwbPh+LPh5KfD6W/fTD5OfIweBbzc1LVA5PNO0OuQjyV8imRSNT9k827Qz0LG90UfaL8feZp2BgZEuXgVMkmx2XT45Oy2fEh2fr4umxofBZM0zIwqcA0C6Z1MA2N37Emd5lhGJHOOvM8ddKOQjfN3aqaehZNwCZlrMTtBofrMMzp1AbpEAs8ZVydADAtmpoqXMAi7ThPFVdxqxSxep3ONFf9kHKK1et0uMtPdq9HknE2QqMhhb2ZGWf52dfH6vVI09wprg02LxfkcPw0iGbZ+zIO/akUsfjPBQ7zU8nmdQzt7W/VTWFW0lT1llrW1vz86Z+aQrK56letSJq26qGxR53E71ppalvzc+xRf+cwuNpK6LV1vtLFPDd6h+dDiduWeLa5H+okz9bitp3jSod5tpJiPdQP9U5+Dg4+Gqz+RmFTi2cr0o3d/ExeJavb0Wa1NtEeys+r1cPkWP/go0fgWEDyEDb9yMejsO2/qv3p1vKzrfOy2jwHV5M0pdWxsf7R5NVY/89jybErYJ4Efx6Ojf1ch2ev5Oej5Cq8/+jwanA0uTrWP3aYPFwFmsn+sdGrw6vD0ZF6PFvKTxb9maQzcPVncN3q4aPk74Bn8tHV4aNR8OpYf/247ZH8HDwcQ1kIDjwcRTxXCzyT4M/6PHunvx0bS16t9v+8Ora6Ogo8V5PJJIQxHbVjyeTUPclPaOoodLGD/auHg4fg3cP+w8PDwUFIUzjsH4X92h/uIT00SFNhxALsD9LHcDjYn9uvybNH8rNFPJj5Su/kZ2s8m45bqU5X+EgPzVea5onO8uhBng8lbh9KP/RQxpVu0gmDnUMXrcdfjTaF1eaqj3WNjp8c10knJ3SRSafOORnRTUxKdeOTtEkKpgiYnGBCZlRzWiSTNlaT+ZuttHR8/3/axrJZqKdZ+6qJa6H9H+McXeckYet+47J3FqFQaIz+ygVRZyQy64w4O3/agPxXj5CB5VsWvu4mIvR5bx2/5EE9/p1FWMDWb6wT1dE8Ox1Jr/64JSyFcTHC9sUHQyyc96Y7MAlvYOsrlu9p1Pnz+6Tv/DdZIuw+7uSXVkACNDv6045XOpOB/VdWs1QkknXu5CXJq3/Yq7NEuD5lMUsjneT5zbXxdprg0nc1Pitv7+mVSlHHBJHu/yBkLaXEjCiGTSXco6JbZzBtvjZCvd6hc/Nfvdl1WITG5VKelkXEzVM8jpqyr275/Gx7eXZK90mXjfY/7gpNDqHRz+So3W60ZO12ockktPjBzXa/KRozwhBTpT8yyIZEk9O1Jw2TQ01ho7nqjXXOk38CZeBfy24ZHcLdxWwUaC7uZg8svx1kPR6PxRGNWfyOxeziPvwEpmzlvX+UEenIeu76i49u4LGS8bSWL24GTda2NbL08eYA5aA/ZnGYYn6HULhvEZqWhcaYPSa0xDye62uL5yALroxmaU/7v1dWBikTtxLn+c0ncy0tvKX7aa2C3zkobPXT+dX3jMwDhte7MT9k6DO70AIvWfsB8PV4FncXFz1AX7jF8BSaFr+pGGLoc1gN/7tS7Tlr9BN2O8qTb6s30Zb/lh9MgKdx7TfLmtG0ZgJnmuxr9jWTZ9fj2Vo0eqK7fuPyYqGXsr+7+Xfp+3L8VPGYNZr30kt5OU8FoJHWKwqbehVtdVSM4V1BGdh3aa6etRgyeWIxjyW2HLN4PKbd2LLJ4nB4IKoLQ87im1ykSKXFkKl8bt42w/akjKfL6/W6FC4Xc3B764kMfCZsK34wXLFT4FkzPw0fRUtGRyNTcgMmejHC5tqT34dSUtsSY/40SO51OZOv6wy3hW3zwMLe9srA0vkCw3PgpbKEp8KmD7u0NiKsUIgVxAY4jFDARgwbAlUSMzXhPT2qEGbMsBs+hS39lldRqNaAP185asi8vJcdtwheGFkjhr7c1GIqAnIowpDae7p0vvL23Hy2tG3eO8s9I1BXxvNEoQh7w66wVp9xiWwuvT4j1mq9Ln5GryXCeq8NeUyvtxGnJ3oX+BPMYUKrP8lM2Wx6LwExILKFvXpbiUNr5afUcxuDUtSQgpbvJid1+atup2cj9LNXB56eLzzdO1sZeIqeKZzzp/kmT75Y77K5bHwXIePbCIU+fEqItS6twgU7Ci+MKl4w2sB46sqEtYqw1mZTADdigyAg0MHRevShEn/W6G9/XWyEZyW7gja0fzs+Pq7K8VRNqHI8F1bMK8BzuyZP4tRlI8B1wNMFvgxDiHpPIBUVU1rtKcEnNsSKTEYLRpvNtqHVe7VhBR/i1ubVgj8JL0SuwuYq4VkjP9WGWP24rUA0Frumd6ArUufidmpqwgBd6t+ZGN1e2qZ5Lj3dM98atzatyxYWg7tOCVkYXAM8T20ZRTh8CnkITOBFa9Pz+V5XxpVRiMMZ1HF5xWGFDXKT8Cq8fIW21J+1+9tXi6h3sYDU2wIv+f1CiwkJPT/Tr/rtFpB7FiO8Z0J2i2XLJFw2CR3Io/6vGF00JJodz32JKkfq3DywtzRgBomwZ2ZMfzWU8UQOIiBuTzMywps5yZycuvQK4pQ4zXjDGf2JHiWcNqMnNrRaWzjDBzPEc+YUivZEj/pf8DJUEzeWnwjOZ8BtP7u2uxwTLmaXo9Go0OM/yGaRy66Xs/uWtbXo8oHD7we75xpGUPRD7MOPcyDN/WWDsqAZlAtVVAI9rkjKxhUxQUAPKoa+1IUO+C4x2vKh2S4Ceig6HOE9hYKAPjkMgwsyi11iBaobprtkAr3f8LhC++NLkDnXUeEykjxr9qwwZt+3+59BHq4JjQ67Qxi9FtrX/NdCT3TRgyoaQeVuiarOW1TV5BDohJ/67qqHoF+2NfDBBnj2GWCecu0XZj0xxHPZnrU8290F5YOk3zJI3EXguc/wtMNGGD0QZoeqp4P838zVvPmp8s48+XyiMT3UwD9YDN98mwWefofRsm/y73uEIP0OQBzE7CD9YsKtXaFn2b9r3PUsRoXZLZiF/nj7yqrypXmlHOZ/zc1X7ubPwqY2GuGpi/R9vSjcNXocax6hHSI2CkmJehrH7hpy6m7MYTfFQPctrjlgirovqjWlVupkExPTsxOT0+sTQ9NDE3pdrgXN8AwTeZoZ1NfYalTN1WtoXgYY+t5ycwAxXZssZcsLwgMYhq7veufZOjyRss8VhULryh+GxbDVi+vp/oZ59vX9UKEaoo618rXca5P9u9vWTVrkSXhB7Wm9epB+oARGvK5TkA5eLf9E7DoFMej1nohrfbyhuM1B+qzqCnUJjNctLE7V5nniEmeQtLOdnIhPwqAJT4iMC6zaMGgDPRKERG1/NrG6LJV+WZOm/9mdnVmPJ2gEPj/jBe1+4gIdC1rJxdeKFS5QTSB/tGLtqa2mP5uIWwT5v2zdytLyZWv3hK7JU+wFUuDTsIvmqQ8jngS4l+YJkxZCX2Oe2ti6SSkMzmxFf8Rg66sWV+Nrx61Lm4EQtXnDiKcro7eBnNVnvISWyJxoT8KnNv3NqXU5mshPBurH1Vxqyrb88OPaPBVhGxQQdTCkEITYRRDokABJD1M1eLHVlgtN5WcOXzsqXBr9U+ur4/XGTwW/QhAoqu5W53mHU7oMP0TLWNq/a8fd2ju8rnmn/+a9WitOTI1RZ1v+T8bxumZ1GCbzq532tTb9Q4fzderqkH+FVIPxP8bb9X8rrUIMQyBdci+lhV9h5jdcE5W75GeO6I8HRv+79p3zMa1//JFeC+Xk8Yn+o8daKHr9Y1QYc5np8anox8Zq5v7mRCtNc7bzqRgGg1wORWKQ0CV3yBRJzlw0RUQ6eWM1mb/ZUtNa+pFaBIu3q+mb4urpmgblumhdyda3Szl7xPy0iM0bvslHaq4cdBDMfe1YunugZIPFW66VwzDC4v3eDPBVG+x8VQXYvK8dOpew5VnKHYFOTGXrzs1S+p6s7HxXBTp63lsZpCIueU6Lptjqhuh7XG5w1OMaRCqWvmmInXO3b4FBxNZzEuUStcog5+rOdnI2e8C23j9B0hxmDU1Vbyns2nod3b//c1MINFc91TX3w9DMYJ3DzMetJFhbr6PTYILOAQt0zXWRZTzBBdVaW+OoHs+uuc61hCeGk5SvSmt9ZdSq1RDgeHWeH3fNdcslPPG4JpCgkE8LfkVJhmsw+pjZ4gkBlnM8VtiZCZC5fawsJLBuzE8sRM1geFrgo3A8QIKLMBz3kQEfeNBHUTjme4/0gc8T9H6AEuDMDgkuptLAk973+QIlDm8tbjuVn3EcPIP5XoR8cSoUClAzITKdIlO+uO8FFQhhcTIFvsRFVCgNJg0VCpApQYIKBHwa8gWJYWBMUek4mSjGeYv52SGeCRwLxRN4YoZKz+AJ4Bkg0/hMKJTyJWYw4DXjiyN/CmZS4OhQIJAiBVTKRyU0JIpbDD4g0KTQr1Pk2Y1xK0gjR2iAZ6AKz3iohCcVRw4m05pQiiR9qTxPDPHEynh2YT+E+RKkj4rjL2bwOA4ZB9FJpkN4nIzTPPEEGaLjFrKRCuCaUJrCE74Eys0QngKepAYPBdLlPLtzXPGF4gHogzCM0oRwQUBD4emQBoITTAIKI0MBDfS3gXQax0MhisTTKR9Tk0oH0LhCpgJQFyOLY0yX6oTcWJHTAtCDQhTndQEML3iIwgo6AaNrF2sKqimI7szPKiBL2kyFqFpVqwH7uGvun1CNZ8FlZS7CBFhFlXo8uzI/C63LjfS+W8ngOfVXVQ6X8eyeedkNHY8Qwmm9l84d5jKRVnVMrZAvTVehfFhR7BVkYMnf69L8JMkATgXAWSSj90iKAsmHspQEjQ/9Lj1kUFQaB4UH0i+eRi8YqEN4AwYbjAyQ5Ty705+JdCgBwwmMnGkqhacELwLpBAXDJy5I+URUXEPBHgZjygswpXxpMgUSEdeQoTgMuwENmYKhVlMa7d2anwkcj8/4kH4VzMRh3E9heGiGYniCTArMgGQCxSDQgCkBrguRSMMHQuQMcAbdDz9PKF02X+nO/hYRYnhiMwkKeM4Az/cYnvEZKs8TS4MJpxKBkA9Jfpon+B/DUmmKIst4dmd+Fv1ZwjMAXAIvijzTFEnHLSijAECgCdE80xDDoTSjJYpx25V6SAAaLoShnigEaegL0Qeg40gQ8TjodhKjQNPhaVC9IQEVB/2nwTUgdkHsIRkYJwWBeKB8XOnO/MTy2i+/OsCog/zwktd3WNGYe0uA5QeeG+NKF85XOoBuzc+28+zKdZP7zbOT6/FdFLeJ95pCqLnqmq7ph5xOpW4iopRCiUxIoeiUzgkw5c15k06JTCqRStdYTeZvtuQR7p9zyg44e84py8+35/A5bZINFnly9pxTlp9vz2F+ymVs+pO75yg+lPyUs3lCD3fPrX0w+flQ+tsHk58qNvOTOz2knmXxvrtcPod4RCpRGuQGpYQuSqVcDoeMKW+WI7OkYL5zTUmb7z3YFM8N0fSITCUbmVZtQBkZUalGNlSMKW9WIfN0wXz3mu08+fb/AZJr89OcuOZgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(filename='data/index.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add required packages. \n",
    "import numpy as np \n",
    "import pandas as pd \n",
    "import matplotlib.pyplot as plt \n",
    "import seaborn as sb #visualization package\n",
    "from scipy.io import loadmat #load .mat files\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_centroids(X, k):\n",
    "    m, n = X.shape\n",
    "    centroids = np.zeros((k, n))\n",
    "    idx = np.random.randint(0, m, k)\n",
    "    \n",
    "    for i in range(k):\n",
    "        centroids[i,:] = X[idx[i],:] #Chooses init centroids from data randomly\n",
    "    \n",
    "    return centroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the function to find the closest centroids. The input to this function is the dataset X and the centroids. Write a line of  code in order to **compute the distance** between the data point and each centroid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_closest_centroids(X, centroids):\n",
    "    m = X.shape[0]\n",
    "    k = centroids.shape[0]\n",
    "    idx = np.zeros(m)\n",
    "    \n",
    "    for i in range(m):\n",
    "        min_dist = 1000000\n",
    "        for j in range(k):\n",
    "            dist = np.sqrt(np.sum((X[i,:]-centroids[j,:])**2))\n",
    "            if dist < min_dist:\n",
    "                min_dist = dist\n",
    "                idx[i] = j\n",
    "    \n",
    "    return idx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the data and call the first function to form intial clusters and visualize the intial clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0xa82a648c>"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAD8CAYAAAD61pSfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XtwHVedJ/Dv0cNWFByNpQhILD+U0UAFghVsxYk3JDCZqYGxB2uq0JikigUCSzyUHVjAs0CqJlub3WUyize7kJiKXZPJkGU2rMsz4AxxssNziFnHRjLIkGSIDbJl2TyUq8Q2kW29zv7Rt6W+53b3Pf263X36+6lKKffevt3nSsn99e93fn1aSClBRERE/hrSHgAREVEeMGASERFpYMAkIiLSwIBJRESkgQGTiIhIAwMmERGRBgZMIiIiDQyYREREGhgwiYiINDSlPYAsufLKK+WqVavSHgYRUa4MDQ29JKXsTHscSWPAdFi1ahUGBwfTHgYRUa4IIU6mPYZ6ML4kK4T4HSHEXiHEvwohXhBCrE97TERElD9FyDC/AOBpKeWAEGIRgNa0B0RERPljdMAUQrQBuBXABwFASjkFYCrNMRERUT4ZHTABdAMYB/CoEKIXwBCAj0spX7U3EELcBeAuAFixYkUqgyQic01PT2NsbAwXL15MeyiRtbS0oKurC83NzWkPJRXC5PthCiH6ADwL4GYp5SEhxBcAnJNS/qXb9n19fZJNP0QUp5GRESxZsgQdHR0QQqQ9nNCklCiVSjh//jy6u7srXhNCDEkp+1IaWt2Y3vQzBmBMSnmo/HgvgDUpjoeICubixYu5D5YAIIRAR0eHEZlyWEYHTCnlrwCcEkK8sfzUHwB4PsUhEVEB5T1Y2kz5HGGZPocJAHcD+Ptyh+wvANyZ8niIiCiHjM4wAUBK+WMpZZ+UcrWU8k+llC+nPSYionp7+umn8cY3vhE9PT24//77q16/dOkS3vve96Knpwc33ngjTpw4Uf9BZpzxAZMo1yZGgEc3APd1WD8nRtIeESVobk7i6z86jXc/eAB9/+WbePeDB/D1H53G3Fy05szZ2Vls3boVTz31FJ5//nk8/vjjeP75ytmpRx55BEuXLsXx48fxiU98Ap/+9KcjHdNEDJgOu9aP8YuJsmXfVuDkD4C5Gevnvq1pj4gSMjcn8edfGcI9X/sJfnL6LF767RR+cvosPvuPP8Gff2UoUtA8fPgwenp6cM0112DRokW4/fbbsW/fvopt9u3bhw984AMAgIGBAXz729+GyVdRhMGA6bC240KwLyae/VPSTh3yf0zGeGL4DA4cfwmTU7MVz1+YnsUzx17CPx09E3rfp0+fxvLly+cfd3V14fTp057bNDU1oa2tDaVSKfQxTcSA6UXni4ln/5S05Tf6PyZjPHJgpCpY2i5Mz+JvnuEJedoYML3ofDHx7J+S1r8TWHkz0NBk/ezfmfaIKCG/PHsh0ut+li1bhlOnTs0/Hhsbw7Jlyzy3mZmZwdmzZ9HR0RH6mCZiwHQYKl0W7IvJ5LN/lpuzob0buHM/cG/J+tneXfs9lEtXtV0W6XU/N9xwA44dO4aRkRFMTU3hq1/9KjZt2lSxzaZNm/DlL38ZALB3717cdttthb/uUsWA6bDlYFewLyaTz/5Zbq6/MCcpPLExxoff1o3LmhtdX7usuRH/7pbwJ0tNTU146KGH8M53vhPXXnstNm/ejDe/+c2499578cQTT1jH//CHUSqV0NPTgwceeMD10pOiM3ot2aBCryU7MWIFlFOHrCyzf2c+MgG/cd/XYQVLW0OTdTJByXl0g3VyYlt5s3XyFvd7qK5eeOEFXHvttTW3s7tknzn2Ei5ML8xlXtbciFt+70o8/L61aGhIP+Nz+zxcS7bIgp615zUb8xu3yeVmN1nI1MLMiXMe3RgNDQIPv28t7n/PW/CWZW248jWL8JZlbbj/PW/JTLAsuiIsjRecHUiAhUDid9ae1y8tv3H376zOPk0W9G+ehOU3VmaLOicpYd6TJ3mt3oTU0CDQf/0y9F+/rPbGVHfMMN0EDYB5zcb8xl20ZpMsnPSEmRM3eR4dyG/1hozEDNNN0LP2vGZjcY87TDaQlQwiC5mafZKS9HvyJAsnMkRlzDDdBD1rz2s2Fve4w2QDWckgTM/U8iqv1RsyEjNMN6aftSclz00r/JtnU16rN2QkZpimqUe3p9cxwmQDzCDIT16rNxnzoQ99CK997Wtx3XXXub4upcTHPvYx9PT0YPXq1Thy5EidR5gPDJimqUeJ0+sYbFrJjixcJkPBzc0BR/cAu94OfL7H+nl0j/V8BB/84Afx9NNPe77+1FNP4dixYzh27Bh2796Nj370o5GOZyqWZFX1bEJJ4lj1KHF6HcOrrOn3OVkKTUYWLpOhYObmgP/zPuAX3wWmJ63nXh0H/unjwPP7gM3/C2gIl+PceuutvjeE3rdvH97//vdDCIGbbroJr7zyCn75y1/iqquuCnU8UzHDVNWzCSWuYzmziebWyteSKHHqllHtcX3x+mw09hRJVuaGSd9P91YGS9v0JPDz7wA//YfEDq1z+y9iwKymfrGc/EF1WSuucpfbscLsyxl4L50DFl+RbIlTt4zqzHKckvzyjvq3Sfv9ceHccP4c3FkdLG3Tk8DBh+o7HqrCgKly+2JRM6O4MkO3Y4XZlxqApifja5JwCwC6jRhegTHJL2+3v02QIBb1b8vLZCisczUyulqvR6Bz+y9iwKzm/KKBsnajHQCClLv8vqz7d3ofI4gkswm/AFArEFWNQyT/5e2Wte/9cOVn2HWrd9CMWsrMSimU3aX5c0WNAFXr9Qg2bdqExx57DFJKPPvss2hra+P8pQsGTJXzi2blv6l8zQ4AQQKUX8Bp7/Y+RhBJZhN+AaBWNqWO62M/Sv7L2+33d0Zpkb90zjvzi3rywVIohbV+a3UPgq25FVi/LfSu77jjDqxfvx4/+9nP0NXVhUceeQQPP/wwHn74YQDAhg0bcM0116Cnpwcf+chH8KUvfSn0sUzG23s5VN3ey9nd+fpe67lfDVf+e63u1lq3yUpqObm4OnD9bh+VxVuATYwAX3wrgBr/XXuNNervLStL/VFm6N7ey7VLFrCC5e/eFqlLNk5Fvr0XLyvx47zkwRk4zgxZgUMnONRaozTMZRU6lwzEdVmB30orWVh/VWVn7c5xXb0WKB2zMkubc6y1glyQIMjLZCishgbgvV+xumEPPmTNWV6xzMosr3tPJoJl0fEvoCvs3FTQcqlOg4rOWOKaS/ObC6tHY0mYrlN1XAOPAFu+7z3WWqXlrDTypCkr3b+ma2gAVv8ZsOVfgL84bv1c/WcMlhlhfIYphDgB4DyAWQAzocsGYbOpoBmHmhnuutUqzzgzG52x1CP7C5tNBcnYwmTKXuPyel+tk4usNPKkiQshRCKlhBD5vwF00afwinLa8vtSyusj1djtrEU0Wtc5njrkfqbtPBPffZv1T5CzcvXL+NK5cEvQqdu847PJZQhBs48gGZtb12vc46/VqMNGHp40RNDS0oJSqZT7YCOlRKlUQktLS9pDSY3xTT/lDLNPSvlSrW2rmn7c+DXBuL3upG5rc2Zcza2Vc21OUZpqqsYlrLm+OJpSav1OVF7NQm6Zp9fiB7WOEUScc5hhts8Dnb+xiZ87BtPT0xgbG8PFixfTHkpkLS0t6OrqQnNzc8XzRWn6KULAHAHwMqy2yV1Syt1e22oFzFqdoerrTl4BT/0yWnyFVYZVg2eUIOE1rjgCT9BuWa8vX7ffw+3/G/jeX1UHzSx05HoxMbjojDfoiRMZoygBswgl2bdJKdcA+GMAW4UQtzpfFELcJYQYFEIMjo+P195blHKd2plplzFP/r/K7eyVevwaVYLyGlccpbWga8uOPgssek35SQFMX7RecytHf++vrC/dlTfrHUM9VhpNKjrlyyBl6Sw03OgshMCyLRnO+IAppTxd/vkbAF8DsE55fbeUsk9K2dfZ2em/s4kR68sdACCsyxXUIOacO7x6rfVPrc5M9ZpBOxjEuVqLPS5VHPNxQdeWlbPA1G/LT0rrMp19W93HYn/pBu3ITbOzVecEIkhwyUuXLud6yXBGl2SFEJcDaJBSni//+zcB3CeldL0xnGdJ1i5HqWVBu3QapqTmViJtaKreV9ylu4kRa6m4M0MLz1291rr0IumSYK1y9bZBqys4jjJ0mosqxF2+zOICEW7yVmam2LAka4bXATgghBgGcBjAk17B0pdX44lbByugV0JTz77thRDUTDLu7KK9G2hWutzsDC8Jfrcec1p+ozW2uMrQ9cp2wi5Or1Yipi9W7qPet2yLA9evJcMZnWEG5Zlh+mVGNudZf5xNH7qZaBBe+0wia1F/F6LRKskuPBG9W9ftdwkkm+14VR3CZMRu/70A7o1gzNwog4qSYRq/cEEs1EUAIIDFS7yXWtOZn9K96L/q2KjMNsOUK9326Za1xFFic/vsK2+OHsi8LsVx/l6S7NCM816fOv+92I1glDyWlsmD6SXZeLjddcOvdBhHOdDZUWrfEDqOW4EB1livXlt+4NG8BMRTDlY/+4qbapftdEra6k2znerRnRnnvT7d/nthA0168tJkRXXHDNNhWeu09QWt3p3E6yzTK4PxW7BclzODuXTOvUwX9ku0vRu46zu1t4vjMoEwvwudZdj8xlKP4OJWdbBLy0F5/Y6i/jcUVtEzLF4eQx44h+kwdPfr5dqOC+4v1vsibLfOyG2D9f0iS+tCdJ2uUK/FHur1BW9yUCn6AgRF//whFGUOkwHTYebeNtnkVaSO2hRTcTmHAK5e438pRxb+p7WDwuizwKLL6xeQ8rhSTtbGE0VeLmNJikl/yzphwCygRDNMtzVm/fbp9T9tGv8z1zt45/ELKwsnOHEx6bNQXTBgFlD/rb1y353L9Ocw/ahf+qcOxXMpR5gvs6gBKOmMI8sBMuzlP3nOyrL896BMYsAsIK3F13W5zbFV3YVEuQZR54sqyJ0+gHiuFUw648hyRqM7tix/BqKEFSVg8rKSpKiddVOvOi7lsMnKtnWddnavyw3c3hvXtYJB13ENKstdibpjS/p3RESp42UlSVEvO7CvPwSqs0T7S1jny9nrEoQgQSfoZRe6iyyEpf6usnTNoe7Ykv4dEVHqmGEmxS/j8MoSdS5W91qvU+fid4hsZj9Zzs6yPDYiqivOYTrEOofpJ4kO2DTWUyUiQnHmMBkwHeoWMLOK3ZFEFEJRAiZLsrRAbRzadav7Oq5ERAXEgEkL1EahS+e48DQRURkDJi1wazLK0iUeREQpYsCkBf07rQUWnLJ0iQcRUYoYMGlBe7f/fT6JiAqMCxdQJV6AT0TkihkmERGRBgZMIiIiDQyYREREGhgwiYiINDBgElH+TIxY9yC9r8P6yRWpqA4YMIkof3TuHUsUMwZMIsqfLN90nIxlfMAUQjQKIX4khPhG2mNJHMtUVBQ6944lipnxARPAxwG8kPYg6oJlKioK3tibUmD0Sj9CiC4AGwH8VwCfTHk4yWOZioqCK1JRCkzPMP8ngP8AYM5rAyHEXUKIQSHE4Pj4eP1GlgSWqYiIEmNswBRC/AmA30gph/y2k1LullL2SSn7Ojs76zS6hLBMRUSUGJNLsjcD2CSE2ACgBcAVQoivSCnfl/K4ksMyFRFRYozNMKWUn5VSdkkpVwG4HcB3jA6WRESUKGMDJhERUZxMLsnOk1J+D8D3Uh4GERHlGDNMIiIiDQyYREREGhgwiYiINDBgEhERaWDAJCIi0sCASUREpIEBk4iISAMDJhERkQYGTCIiIg0MmERERBoKsTQeERGFN1qaxPa9wzhy8mWsWbkUOwZ6saKjNe1h1R0zTCIi8rV97zAOj0xgZk7i8MgENj74DEZLk2kPq+4YMImIyNeRky9XPD5/cQbb9w6nNJr0MGASEREAq/S6eddB9NyzH5t3HZzPItesXFq1rRpEi4ABk4jIEF4BT5daet2+dxijpUlcnJ6t2tYtiJqOTT9EVCh2A8vQiQm0Lm7C5NQs1kZsZEmrKUY97qXpWQyPnQWA+YC3Z8t67f2pWeORky9j+95hHC3vEwAEgBu627FjoDeWz5AnzDCJKH0TI8CjG4D7OqyfEyOJHcrOomalNRc368imou5zJoZ9OdXKGNXjOgMbELxsqmaNa1YurdpHY4PAni3r2SVLRJSKfVuBkz8A5masn/u2JnYoryASJLiogWzIJTMLs59aAfHtn/9uxXbqcaSy/6Bl0x0DvVjX3Y6mBoF15SzSLYgWFQMmEaXv1CH/xzHy+sIPEgjUQDY7VxmqdPdVKzN1C4jO7dTjrO5qqwp4QazoaMWeLetx/HMb5rNItyBaVJzDJKL0Lb/RyiydjxOyY6AXGx98Bucvzsw/t6SlKVAg8Msge7vatPflNmfotGblUhwemfB8346B3lBzp0HmXO0gSgyYRJQF/TutMuypQ1aw7N8ZaXd+AWFFRysuTFV2fV6Ymg00J+cVyABgcXOj9r7U/agZox0Q1WPZ24UNZs59Hh6ZwLbHj6ClubHwK/nUwpIsEaWvvRu4cz9wb8n62d4daXe1Sp1R5+XsMqWbIHOhXuVOe27ztv/+PQDAVz9yU9V2US4hUcd4dOxsoKYl9dhNHcsXaR88x5hhEpFRRkuT+KGSkakBwq2UqbNft/eo5d0gwdcrQ1QzwAe+9WLVdpt3HazYxu0SEq9M2y9DBmoHfXV8V278RLQznJxghklERtm+d9i1W9SZFW3fO4wdA70VzS1+RkuT2PjgM1VZ2IqOVjx59y2xN8XUmtvU3cYr0/bLkIHaQV891qLX/e7lvm8wBAMmEeWWW1nSLXA45wLDXCu5fe9wRRYJWJnV5l0HAaCqszQqnZKx1zbO34lXpm1ntmrQFIBW0FePPfXrn7/q+wZDGB0whRAtQojDQohhIcRzQoj/lPaYiKgshsUK3IKg+mW+rrsdKzpaPTMyv7lA+zWv8mWcixQ46VzK4bWN83eiZtqzc3L+MzqXvBOwunv/5S9+Xyvoq8d+6cn/kdxKExkipFR/peYQQggAl0spfyuEaAZwAMDHpZTPum3f19cnBwcH6zpGosJ6dEPlpSQrb7YafgLouWc/Zua8v8N6u9rw4B1rsKKjtSrwretux54t6z2fB+AbLG1NDQLHP7dh/nHa946s9TsBMJ9Zen3uoIQQQ1LKvlBvzhGjM0xp+W35YXP5H3PPEIjyJMBiBUHuouHkvMTDKyNTM0+71OpV3l3SUtkr6czagOSWydOl03R0eGSi6kRg6MREpIXbi8DogAkAQohGIcSPAfwGwDellMktIUJE+tTFCXwWK6jVvNLUICBc3ucMeG6r2ADuAca+SbKaqa3rbp9v8rGpq+/oNOMkqVZDj5fWxU0Vv+Ntjx9hAFUYHzCllLNSyusBdAFYJ4S4zvm6EOIuIcSgEGJwfHw8nUESFVH/TqsM29Bk/fRZrMAtCxwtTVYEwRtcgoROtuUVYNQmH3sFH/uYTQ2VIdoeo1/DTtTbb+nwauhZ0tLkeWKxrrsdk8piDkGvzSwC4wOmTUr5CoDvAniX8vxuKWWflLKvs7MzncERFZHmYgWjpUlctqix6nn1C3zHQC9Wd7UBWGhi0bnEwyvAqNQVfNTA+Karr7AWYj8xgSUtTWh0adiJo1yrG3TVEvSTd9/iemJhz12uVT6PGliLeMNolelNP50ApqWUrwghLgPwzwD+Wkr5Dbft2fRDlICJkepl7wKs5OPVeKM220TlbNa5bFFjVYYJLFxysaKj1fdelPa2ahON2pAT5jP4NSnpfD71PqCf/MM34IFvvVj1vNfncWtqWnnl5Wz6McBVAL4rhDgK4Iew5jBdgyURJSTirbu8Mpu4bzPlLO8+efctVc09QOVcpTon+tyZc1XbqtmfOma1YUhH2DlS9T6ga1cuxZ4t6/HAt150ff7BO9a4LsfntoBDURgdMKWUR6WUb5VSrpZSXielvC/tMREVjtr9evIHga65VIOM28X1XmXKsHOGzhV8VEECuFvZ2K1hKEiDTdh1cL0Crdfzbk1Sbgs4FKlUa3TAJKIMcOt+dckyR0uT2PTQAaz6zJPo/syT6H/oAEZLk1VzcW4X13vNDUaZM/Sa2/QKUDsGen3n/ZylTHU7tcFGvVG0epwwS/F5BdogAdgtOBbphtJGz2EGxTlMogRMjABffCsqLoFuaLKafRzc5iqXtDThwtRsxWLnbosCeM0NxjFnGGQhAvUzOMevzgk6CbhfIK5+/igLIHh9Dr/PV2uedklLE568+5bCzGEyYDowYBIlRGNVn1or1PitThNmFZ8k6DQOOTU1iJrB1OY3djWw2Y08YVYb8vsMq7vaXO+bWZSVfhgwHRgwiQII0v2qsW2tZejs6x7dMsYw2ZObOJe10zkBsAOg3UzjF2D9smO3zNa5ryAnCn5/B68xFCVg8n6YRCaLeEmHL7v7FVjofvVaC9a+5tLHjoFebHv8CI6OnYUA0LqoEa86Lqa358qcX+b2c173lfR63ot6n0e3e0zq8rvnpPMaUTtI18pG1QUQnIF9SJlbjNKY47dtkeYr3TBgEpksSFALKsBasDpWdLTiiW1vm3/sdcPmoDd+DiLOZe12DPR6Zo3ORRCcQdrJrfxp/06c2x8emajKKNXHQQKdGujd5pGLiiVZB5ZkyTj3dVjXP9pcmm1Ci+FuI/VWq+Qa95znaGlyPmtW2fOXQydfxqzLmrVu5WCvcmmjANauao99DlP3vUUpyTJgOjBgknGSDGpJlnsTogYcNYuLEmiCHNcpyHyj17xo0s1MtTBgFhADJhknh0EtSUEaceLkzNpmlRs7q9mhV5B2awwSAN7S1QYB4Lkz51K5/yZQnIDJOUwik2k02xSJOj+nXv+Y1Ko1zuYjNdtcu0ovSKuNQfY1kHE2KpE/rvRDRIWhrpJj393EVo8u0LAr9ajB/MLULFZ0tKZ+/80iYYZJRIWhXmbi1YlbzzHoUrNj59J2bs9T/DiH6cA5TCJzxLkIQRbEtThDEooyh8mA6cCASWSOei+Ll1VJBlR734d+Pi5P/rdNxk/xGf8BiaiYOLdniXLHFt19i4ZG9QYsRmLAJCIjhb1vpGmSPHEo2kkIAyYRGSlsN6ppkjxxKNpJCOcwHTiHSUSm4RxmfBgwHRgwicyUhU5SkxWlS9b4MwIioiQbX+I2WprE5l0H0XPPfmzedRCjpcm0h0RlDJhEZLw8dczmKbgXDQMmEeVKmAwsTx2zeQruRcOASUS5EiYDy1PHbJ6Ce9FwLVkiypUwGVjY9VvTsGOgt+7r25IeBkwiyhXTFxvPU3AvGpZkiShX8lReJbMYnWEKIZYDeAzA62DdJ3a3lPIL6Y6KiKJgBkZpMTpgApgB8Ckp5REhxBIAQ0KIb0opn097YERElC9Gl2SllL+UUh4p//t5AC8AWJbuqIiIKI+MDphOQohVAN4K4JDy/F1CiEEhxOD4+HgaQyMiohwoRMAUQrwGwD8A+PdSynPO16SUu6WUfVLKvs7OznQGSEREmWd8wBRCNMMKln8vpfzHtMdDRET5ZHTAFEIIAI8AeEFK+UDa4yEiovwyOmACuBnAvwVwmxDix+V/NqQ9KCIiyh+jLyuRUh4AINIeBxER5Z/pGSYREVEsjM4wiUw2WpqsWqR7RUdr2sMiMhYzTKKc4o2GieqLAZMop3ijYaL6YsAkyineaJiovhgwiXIqrttcjZYmsXnXQfTcsx+bdx3EaGky5pESmUFIKdMeQ2b09fXJwcHBtIdBVFebdx2suCHzuu523j6LAhFCDEkp+9IeR9LYJUuUgDx1sHIulEgPS7JECchTByvnQon0MGASJSBPWVtcc6FEpmNJligBa1YurZgXzHLWtqKjlXOWRBqYYRIlIK9ZW62OWXbUUpGxS9aBXbIZMjEC7NsKnDoELL8R6N8JtHfn7xg5U6tjlh215KYoXbLMMCmb9m0FTv4AmJuxfu7bms9jBJR2Bldr7jVPc7NEcWPApGw6dcj/cV6OEZBOd22SQbVWxyw7aqnIGDDzYmIEeHQDcF+H9XNiJO0RJWv5jf6PM3gMv0CmG+R0MrgkL1mpNfea17lZojhwDtMh03OYj26wyoa2lTcDd+6Pvt96zeMFPU4O5zD95vd05/50tuu5Zz9m5hb+v21qEDj+uQ2hx00UVVHmMBkwHTIdMO/rsObabA1NwL2l6PtNKhCndZwU+QUy9TX79TddfQUEgOfOnMOalUvxyT98Ax741ou+KwSpQXV1VxtamhsxdGICrYubMDk1i7U1VhfSWYkoT6sVUbqKEjBZks2LpEqU9ZrHC3OcnJWh/eb33Ob6ZuYkjo6dxfDY2fny6gPfehF7tqzH8c9twJ4t610DlFoWFQAOj0xgVgLnL85gVqNUq1PWDVv6TbtxiSgpDJh50b/Tysoamqyf/Tvj2a8aeJtbkwlMYQJ+BrtY/fjN7zlfEz770Ok6tRcasIPqc2fOBd6X+trQiYmqIBe2IzZPywISBcGVfvKivTuZEmb/TmDXrcCl8pfupXNWYIr7WP07q+cLa6mVlWbsOkp1xRw701JLmmpJ1SlM16m6qpDOvtT3tC5umn9sB7mwqxXx0hMyFTPMOKVRQox6zPZuYFopmbmVS+M4zp37rXnXO/frBbZaWWnGM1CvTMuZba7uakNvV1ukrlN7f40CuHxRIwBAALg0PetZDlWz4cmp2YrXj5x8GTsGerG6q813f27l11qXnrBkS3nFph+HyE0/aTS2xHFMnX2o2yy+Atjy/YXAl0S2V2ufSTVCxSTubladJhy3LtsdA72h3rdny/pQK//UOh5XCzJPUZp+WJKNUz0vhLeDiTOIhT2mTrlU3a9aunWOxc72op4sqGVoO8u1x/n6XuDM0MLrSVyrGUHcC7DbGSuwUDZVA41bOVR937bHj6CluRFHTr4836X709NnsaSlssPWa3+1HtdazF19z+GRCYyWJtmBS5nHkmyconSy+pU83V5zC5ZBj2nTKZe67dcZROtxsqCWYIHQjVD1KAvqlDR1jZYm8UNlntJtbtCtHKpud3Ts7HwFmlG4AAAKlElEQVSp2O7StTts165cWtGdm8TKP27bsDGI8sDogCmE+FshxG+EED+tywGjdLL6zce5veYWkOLsnlX177TKsE7OIFqPlXnUz/yr4eDzomX16ORc0dGKlmZrTlECGB47G/o42/cOQ508cQs8bp26QTJbNbgmsfLPjoHeqk5hNgZRHphekv07AA8BeKwuR4vSyeqXobm9tvzG4HOXUeYZ27utOUuv0m2YLtig1M8cISgH7eS05w+DLA4Q5jhuxzxy8mXMzlX3GjiDkzq/+Z1PvWN+XDsGevH2z3+3KuC6UYNrrfJqmHtpruhoxQ3d7RXl6tk5ic27DnJxBMo0ozNMKeX3Abj372eNX4bm9lqYbDZqV6lf6TZMF2xQIT6zV+lVDQz2F7ZXydTOSNXFATY++IxvmTXKYuXOLFgNduu62ysCi1/GbAcopyUtTYG7dOMsY9uZqU0CvGaTMs/4LlkhxCoA35BSXldr21SXxvPL/uLqQE2iqzTlayGdmZW6zJzdrenWkWm/T71+0atj021pu1rvUcendozW6nr1Wk5PZ1u1IzeOZe6S6G7lurhmYJdsQQgh7gJwFwCsWLEivYH4lXPjWrQgxpLmvCS6YwNwBr2jY2fnn7ezFa+SqF1KVL+wvUqmXosD+L3HeZxaY3frelWP6ReganXkhimdqpJYkCDuTmKiJBldktUhpdwtpeyTUvZ1dnamPZxkJbG8Xsr3lKy1/FtcXZ7OxQEalY6VsF/ytQJQkIaaetx2K+q9MN1KurxdGOUJS7IOmb5bSValfBcSv2Xm7Es6jo6dhSg/fvCONRWlyCClyrCNP7pjz/oF/FHvcJK3z0v6ilKSNTpgCiEeB/AOAFcC+DWA/yilfMRrewbMEOo8h6l+ITtvh6XOYV6ansWwo0wb9Qs67i/8NG6flfQx/X5HnK80V1ECptFzmFLKO9Ieg/GSWhTegzrvZ98Oy03PPZXjijrnFvccXhzzikHprBZUy2hpEtseP+Kaufv9jjhfSXlX+DlMCsfrEoOkV9DxWlbNTdxzbm++unLhhjx+4esGfb+/4/a9w/MNVuqCDH6/c85XUt4xYNZDzm6ErMPrur+kV9AJsqxalC/o0dIkNj74TMVnkUDuv/B1TyL8/o5uQdZ+zu93rt7HkwsUUN4YXZLNjJQvvUiCV6aS9L0Q3Vat8TpGlJLn9r3DOH9xpuK558+cy/Wc22hpEhenrdt42aVUr6AfpLRqPwekU2YmqhdmmPWQ8qUXSfDKVKKWQWtxW7UmidKozsLmSUmqrK2WUhc3N3pmebVKq85F5Xt9Ai+RSRgw66EeC5M71aEE7FV6q8c8VRrXHC5paYp8nNHSJDY9dACrPvMkuj/zJPofOuAaDJMqawfJ/muVVp/Y9jacuH8jRu7fiH3b3sbyKhWC0ZeVBJXYZSX1Xj4u5Wsjw0jjEot6j8ftmlG3S1OSuvyC10FSUnhZCcWnzpde5LEEHMflDnEKMhenG1z9mmWckrr8wl5b1znOpEQ54cjayRORjSVZE9W7BByDpJuFkqRbQnULfLr3tIxDPbtUo5SV63GvUqIwmGGaqB73poxZni9q1w32OwZ6qy74dwuGJnSaRjkByvPJE5mNAdNE9S4Bx6Ce5cK46QZ7u1nGNG4l1CAnQOr733z1FRVLGubp5InMxqYfB64lS2EUfc7NrZnI7QTI63eivn91VxtamhsL+/vMIzb9EBkq7gBnQgk1CrcSapDfifr+vC8QQeZi0w9VSXo92LSxqcRd2L+73yIHOvtMerELorgwYFIV0wNKUZpKggbAsH93v65enX1yUXbKC5ZkqUocASXL83p57sgNIui1rWH/7n7lV519Fr2kTfnBDJOqxFEiy2KWamdcQycmsKSlCY2GZzRBA2ASpdF6lFtNn0Kg7GDApCpxlMiyWPa0g/isBM5fnMHalUtze5upJOYGkyiN1qPcmsWTMzITS7JUJY4SWRbLnlkM4mHplFuDXtuaRGm0HuVWk/6ulG0MmJSILC5EUCuIZ3neVcW5wQVZPDkjM7EkG1UdbqWVR/Vct1RXrfJgnkp7vBRjAbtsqV640o9DqJV+cngrLXKX1G21kpCHbDgPY6R4cKUf0pPDW2mRuzyV9vJQbs3aLduIomJJNqoc3kqL3LG0Fy8245BpmGFGlcNbaZG7PGRteZKnjJ1IBwNmVDm8lRZRPWSxU5ooCgZMIkoEM3YyjfFzmEKIdwkhfiaEOC6E+Eza4yEionwyOmAKIRoB7ATwxwDeBOAOIcSb0h0VERHlkdEBE8A6AMellL+QUk4B+CqA/pTHREREOWR6wFwG4JTj8Vj5OSIiokBMD5g1CSHuEkIMCiEGx8fH0x4OERFllOkB8zSA5Y7HXeXn5kkpd0sp+6SUfZ2dnXUdHBER5YfpAfOHAH5PCNEthFgE4HYAT6Q8JiIiyiGjr8OUUs4IIbYB+L8AGgH8rZTyuZSHRUREOWR0wAQAKeV+AFyKh4iIIuHtvRyEEOMAXvV4+SUAV2rsJs7tsrqvohyT4y/eMfM+/jSO+RKAN0opl2jsL9eMzzCDkFJ2CiFcb4gppezzei2p7bK6r6Ick+Mv3jHzPv40jqm7LxOY3vRDREQUCwZMIiIiDSzJVtsd8rWktsvqvopyTI6/eMfM+/jTOKbuvnKNTT9EREQaWJIlIiLSEKgkK4QQAL4GYCOshQAmAVwAMA3gcgCXlR9fBqA51pESERFFN1f+KRzPTZcfb5ZSft3rjUEzzI0A/gjAGQB3AGgFcA7ALIAWAH8K4Ffl/V4CMFV+36/L28jyYOfgbUpzOzdBt9fBmjURUfa4fTfPKs+PABiGFRBnAPwE1l2rjpa3+0sAfwLgp7CSwH/2O2DQgPlhWIuX/yuA0fIAXikf+CSAXgDfAfDb8msXYS0E8KpjwLL8vE0qH9LOTOdCjC8onWAoam9CRERlcScZc8pPr+PNwQp6s+XHs7AC4SosxJhzAM7DSvbmYMWkNlgJ34SUctJvIEG7ZJcDOAvrHpPLYAXCJbCCypnyc9PlwTWVB3Q5rLuENJcH3KAcV5Q/pPNxmLEBwQMsgyERUbzi/l5tUH6qx1GPZ8eORgDvLv/7TPn5NwAowaqEdgP4z+XtJIBP6Q4kKZcDGIeVkdpZ5QXHcWfK/0wnPA4iIjLbtPLTWblsLD9/FsBVAG6CFZs+C+Cj5e22CCF8Y2LNgCmE2CqE+LEQ4sewyrFtsDLN07AC4vnywa4uP9dcHtzi8vMSVpCcLR+vBQtnAE3lf+zs02aXaadARERUmx1zxsuPneVZOwt9DlY5thHAFQD2wGpSlbBik++6uTUDppRyp5Tyeinl9QD+BlbZ9VoAK2AFu98p72clrMnVm7FQpgWA18IKpq1YCKDO4DjjclhR3uf5WuNT2DVuNuoQEWWX+h2t853t3GYalb0wNgHgNY5thfLvF8r//gqs2PMHALbASvZasBBsXQVauKB8WckTAN4FK0JfKP+zGFa2OQcroi/S3ikREVF9zaFyDvRFAHdLKX27ZLnSDxERkQau9ENERKSBAZOIiEgDAyYREZEGBkwiIiINDJhEREQaGDCJiIg0MGASERFpYMAkIiLS8P8B9w/L5RjUWK4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa88a086c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = loadmat('data/GroupA.mat')\n",
    "X = data['X']\n",
    "k=2\n",
    "initial_centroids = init_centroids(X,k)\n",
    "initial_centroids\n",
    "idx = find_closest_centroids(X, initial_centroids)\n",
    "idx[0:3]\n",
    "ax = sb.swarmplot(x=X[:,0],y=X[:,1],hue=idx)\n",
    "ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.75033133, 4.70376783])"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i in range(2):\n",
    "    ind = np.where(idx==i)\n",
    "    \n",
    "ind1=np.array([0,1])\n",
    "X[ind1]\n",
    "np.average(X[ind1],axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this function we compute new centroids for from the clusters. Write a line of code to **compute the centroid**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_centroids(X, idx, k):\n",
    "    m, n = X.shape\n",
    "    centroids = np.zeros((k, n))\n",
    "# compute the new centroids   \n",
    "    for i in range(k):\n",
    "        indices = np.where(idx == i)\n",
    "        centroids[i,:] = np.average(X[indices],axis=0)\n",
    "    \n",
    "    return centroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the function where we run kmeans and optimize the clusters. Call the function we defined above in proper sequence following the flow chart above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_k_means(X, initial_centroids, max_iters):\n",
    "    m, n = X.shape\n",
    "    k = initial_centroids.shape[0]\n",
    "    idx = np.zeros(m)\n",
    "    # 1. Call function to initial_centroids\n",
    "    centroids = init_centroids(X,k)\n",
    "    for i in range(max_iters):\n",
    "        # 2. Form intial cluster by finding the distance from centroids\n",
    "        idx = find_closest_centroids(X, centroids)\n",
    "        # 3. Update the centroids\n",
    "        centroids = compute_centroids(X, idx, k)\n",
    "    \n",
    "    return idx, centroids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main function where give intput for no of iterations for optimization.\n",
    "initial_centroids = initial_centroids = np.array([[3, 3], [6, 2], [8, 5]])\n",
    "idx, centroids = run_k_means(X, initial_centroids, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xa796290c>"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAHVCAYAAAATqShMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X+UXGWd7/vP0+kiCaSNGQgehHISnBnILzpDmm4xF8jF648AhyPGpRyJQ+IAySiN3NFV6xxdWffeybq66OsacSJHAqMnQuJRVxzENZIZZo4YHBZW0o2ddEiAo5GhCEISljQdgdCdeu4fO9U/qndV7arav+v9Wiur6eqqXU89tZv+1nd/n+9jrLUCAAAA0qot6gEAAAAAQSLgBQAAQKoR8AIAACDVCHgBAACQagS8AAAASDUCXgAAAKQaAS8AAABSjYAXAAAAqUbACwAAgFRrD+Kg55xzjl2wYEEQhwYAAAAkSQMDA8ettfNr3S+QgHfBggXq7+8P4tAAAACAJMkY8+9e7ueppMEY805jzE5jzDPGmEPGmMubGx4AAAAQDq8Z3m9I+idr7ceNMWdIOjPAMQEAAAC+qRnwGmPmSrpS0jpJsta+LentYIcFAAAA+MNLhnehpGOS/rsxplPSgKTPW2v/MPlOxpjbJN0mSe95z3v8HicAAEBijY6O6sUXX9Rbb70V9VASadasWbrggguUyWQaeryx1la/gzFdkn4paaW1Nm+M+Yak1621myo9pqury7JoDQAAwPHb3/5WHR0dOvvss2WMiXo4iWKt1auvvqqRkREtXLhwys+MMQPW2q5ax/CyaO1FSS9aa/Onv98p6dK6RwsAANCi3nrrLYLdBhljdPbZZzeVHa8Z8FprX5ZUMMZcdPqmD0g62PAzAgAAtCCC3cY1O3deuzT0StpxukPDYUnrm3pWAAAAICSe+vBaawettV3W2kustR+11v4+6IEBAADAPy+//LJuvPFGvfe979XixYt1zTXX6LnnntPzzz+vpUuXNnTMbdu26aWXXmpqXM8884wuv/xyzZw5U1/72teaOlYlngJeAAAAhKhQkHp7pe5u52uh0NThrLW64YYbtGrVKv3mN7/RwYMH9ZWvfEWvvPJKU8dtJOAdGxub8v0f/dEf6e/+7u/0xS9+samxVEPACwAAECeFgtTZKW3dKu3d63zt7Gwq6H3ssceUyWS0cePG8duWL1+uK664Ysr9tm3bpttvv338++uuu04///nPderUKa1bt05Lly7VsmXL9PWvf107d+5Uf3+/brrpJi1fvlxvvvmmBgYGdNVVV2nFihX68Ic/rN/97neSpFWrVulLX/qSrrrqKn3jG9+Y8pznnnuuLrvssoZbjnnhtYYXAAAAYejrk06ckEZHne9HR53v+/qkLVsaOuSBAwe0YsWKhoc0ODioI0eO6MCBA5Kk1157Te985zv1zW9+U1/72tfU1dWl0dFR9fb26uGHH9b8+fP1gx/8QF/+8pf1ne98Z/wxu3fvbngMzSDgBQAAiJN8fiLYLRkdlfbsiWY8ki688EIdPnxYvb29uvbaa/WhD31o2n2effZZHThwQB/84AclSadOndJ55503/vNPfvKToY23HCUNAAAAcdLTI5Vf3s9knHreBi1ZskQDAwM179fe3q5isTj+fan37bx587Rv3z6tWrVK99xzj2655ZZpj7XWasmSJRocHNTg4KCGhob06KOPjv/8rLPOanj8zSLgBQAAiJNcTpozZyLozWSc73O5hg959dVX6+TJk7rvvvvGb9u7d++0EoMFCxZocHBQxWJRhUJBe05nlY8fP65isag1a9Zo8+bNeuqppyRJHR0dGhkZkSRddNFFOnbsmJ588klJznbKTz/9dMNj9hMlDQCSoVBw6tfyeSf7kctJ2WzUowIA/2Wz0r59zv/z9uxxMrtN/j/PGKOHHnpId955p+666y7NmjVLCxYs0N133z3lfitXrtTChQu1bNkyLV26VJde6myue+TIEa1fv348+/vVr35VkrRu3Tpt3LhRs2fP1pNPPqmdO3fqjjvu0PDwsMbGxnTnnXdqyZIlVcf28ssvq6urS6+//rra2tp099136+DBg3rHO97R8Oud9vqttb4drKSrq8v29/f7flwALaq0Yrm0iKOU7di3j6AXQCIcOnRIixYtinoYieY2h8aYAWttV63HUtIAIP6qrVgGAKCGdAS8PjdnBhAzMVyxDABIjuTX8JZf6hwclHbs4FInkCY9Pc7v9uSgt8kVywCA1pH8DC+XOoNB1hxxEsCKZQBA60h+hpdLnf4ja464CWDFMgCgdSQ/4OVSp/8C2NIQaFo2y/kHAGhI8ksauNTpP7LmAACkzssvv6wbb7xR733ve7V48WJdc801eu655/T8889r6dKlDR1z27Zteumll5oa144dO3TJJZdo2bJlev/73699+/Y1dTw3yQ94S5c6N2xwsrobNnDpvVkBbGkIAAC8KwwX1PtIr7rv71bvI70qDDe3lsZaqxtuuEGrVq3Sb37zGx08eFBf+cpX9MorrzR13EYC3rGxsSnfL1y4ULt379bQ0JA2bdqk2267rakxuUl+SYPEpU6/5XJOzW55k3+y5gAABK4wXFDnvZ068fYJjRZHNfjyoHYM7dC+jfuUndtYQu+xxx5TJpPRxo0bx29bvny5JOn5558fv23btm3q7+/XN7/5TUnSddddpy9+8Yu64oor9Jd/+Zfq7++XMUaf+cxnlM1m1d/fr5tuuml8p7WDBw/qr//6r3XixAmdc8452rZtm8477zytWrVK73//+/XEE0/o+uuv1xe+8IXx53z/+98//t/ve9/79OKLLzb0GqtJfoYX/iNr3jy6XAAAGtT3RN94sCtJo8VRnXj7hPqeaLwD1YEDB7RixYqGHz84OKgjR47owIEDGhoa0vr16/Xxj39cXV1d2rFjhwYHB9Xe3q7e3l7t3LlTAwMD+sxnPqMvf/nL48d47bXXtHv37inBbrlvf/vbWr16dcPjrCQdGV74j6x54+hyAQBoQv5IfjzYLRktjmrPkejW0lx44YU6fPiwent7de211+pDH/rQtPs8++yzOnDggD74wQ9Kkk6dOqXzzjtv/Oef/OQnqz7HY489pm9/+9v6t3/7N38HLzK8gP/oDQ0AaELP+T3KtE1dS5Npy6j7/MbX0ixZskQDAwM179fe3q5isTj+/VtvvSVJmjdvnvbt26dVq1bpnnvu0S233DLtsdZaLVmyRIODgxocHNTQ0JAeffTR8Z+fddZZFZ93//79uuWWW/Twww/r7LPPrueleULAC/iNLhf+8FIWQukIgBTKrcxpzhlzxoPeTFtGc86Yo9zKxtfSXH311Tp58qTuu+++8dv27t2r3bt3T7nfggULNDg4qGKxqEKhoD2n/3YdP35cxWJRa9as0ebNm/XUU09Jkjo6OjQyMiJJuuiii3Ts2DE9+eSTkqTR0VE9/fTTNcf2wgsv6GMf+5gefPBB/dmf/VnDr7EaShpQWaHgZCXzeadzA43+HbXmhd7QzfNSFkLpCICUys7Nat/Gfep7ok97juxR9/ndyq3MNbxgTZKMMXrooYd055136q677tKsWbO0YMEC3X333VPut3LlSi1cuFDLli3T0qVLdemll0qSjhw5ovXr149nf7/61a9KktatW6eNGzeOL1rbuXOn7rjjDg0PD2tsbEx33nmnlixZUnVsf/M3f6NXX31Vn/3sZyU5Web+/v6GX6vr67fW+npASerq6rJ+DxQhKw8mSp0aWj2Y8DIvzF3zenulrVunf2jYsGGittzLfQAgJg4dOqRFixZFPYxEc5tDY8yAtbar1mMpaQhDEi+7Uofqzsu8pL3LRRjns5eyEEpHAAAeUdIQtKRediWYcOd1XtLa5SKs89lLWQilIwAAj8jwBi2pmVJ2W3PX6vMS1vnsZctwthUHAHhEwBu0pGZKCSbctfq8hHU+eykLSXvpSAtJYtUXgGShpCFoSb3sWgom+vqcYKa7my4NEvMS5vnspSwkraUjLSSpVV8AkoUuDUFjxT7ShPMZPqPZBloFXRqaR5eGOOOyK9wk9Rou5zN8ltSqLyCJXn75Zd14441673vfq8WLF+uaa67Rc889p+eff15Lly5t6Jjbtm3TSy+91NS4Hn74YV1yySVavny5urq6AtlamJKGMHDZFZPVuoYb9w0/OJ/ho6RWfQFB8/tPgbVWN9xwg26++WZ9//vflyQNDg7qlVdeUbaJA2/btk1Lly7Vu9/9bs+PGRsbU3v7RAj6gQ98QNdff72MMdq/f78+8YlP6Jlnnml4TG7I8AJhq9bpoBQMb90q7d3rfO3sTE4GGKhTq68DBdwE8afgscceUyaT0caNG8dvW758ua644oop99u2bZtuv/328e+vu+46/fznP9epU6e0bt06LV26VMuWLdPXv/517dy5U/39/brpppu0fPlyvfnmmxoYGNBVV12lFStW6MMf/rB+97vfSZJWrVqlL33pS7rqqqv0jW98Y8pzzpkzR8YYSdIf/vCH8f/2EwEvELZq13CT2sYOaBBVMsB0QfwpOHDggFasWNHw4wcHB3XkyBEdOHBAQ0NDWr9+vT7+8Y+rq6tLO3bs0ODgoNrb29Xb26udO3dqYGBAn/nMZ/TlL395/Bivvfaadu/erS984QvTjv/QQw/p4osv1rXXXqvvfOc7DY+zEkoagLBVu4ZLQSNaEFUywFRx/FNw4YUX6vDhw+rt7dW1116rD33oQ9Pu8+yzz+rAgQP64Ac/KEk6deqUzjvvvPGff/KTn6x4/BtuuEE33HCDHn/8cW3atEn/+q//6uv4yfACYat2DbfVN7YAAATyp2DJkiUaGBioeb/29nYVi8Xx79966y1J0rx587Rv3z6tWrVK99xzj2655ZZpj7XWasmSJRocHNTg4KCGhob06KOPjv/8rLPOqvn8V155pQ4fPqzjx497eVmeEfACYSl1ZlizRrr+eulTn5p+DZeCRgBoeUH8Kbj66qt18uRJ3XfffeO37d27V7t3755yvwULFmhwcFDFYlGFQkF7TqeVjx8/rmKxqDVr1mjz5s166qmnJEkdHR0aGRmRJF100UU6duyYnnzySUnS6Oionn766Zpj+/Wvf61Sm9ynnnpKJ0+e1Nlnn934i3VBSQNaV5jdENw6M7j1r231jS0AAIH8KTDG6KGHHtKdd96pu+66S7NmzdKCBQt09913T7nfypUrtXDhQi1btkxLly7VpZdeKkk6cuSI1q9fP579/epXvypJWrdunTZu3KjZs2frySef1M6dO3XHHXdoeHhYY2NjuvPOO7VkyZKqY/vRj36kBx54QJlMRrNnz9YPfvAD3xeusfEEWlPYGyjQXT9acW/1BiD12HiieWw8AdQr7G4IcVyB0Cpo9QYALY+AF60p7ACUxWjRodUbALQ8Al5EL4ptdsMOQFmMFh2y6wBiIogy0lbR7NwR8DYrimAtTaK63Bx2AOrWXX/XLifLyLkTLLLrAGJg1qxZevXVVwl6G2Ct1auvvqpZs2Y1fAwWrTUj7IVPaRTlYq7SQqYouiEUCtKyZdLIiFQsSm1tUkeHNDTEueM3fk8BxMDo6KhefPHF8b62qM+sWbN0wQUXKFOWwPC6aI22ZM2oVhtYK1hj1bij0uXmH/7Q/7lxm/OoOiRs2iQND098Xyw632/aJG3bFsxzBnHOJeE8ptUbgBjIZDJauHBh1MNoWWR4m9Hd7VyGd7s9n6/8uKAyTkkIPsq5ZXglJ+NZLPo7N3HK8r3rXdLRo9NvP/dc6ZVX/H++IF5/3OYUANByaEsWhkZrA4NYNZ7U1kvltbSlRtOlbQ39WlEfxUr9ONV3B/H6Gz1mnOYFANASCHib0ejCpyBWjSe19VL5Yq7586ffx48V9WGv1K/1AWT1avfHVbq9WdVKRxoNPBuZ06R+MAMAJBoBbzPcVt57uZwbxKrxOLReajRzl806tbT5vPSJTwSzoj7slfqVPoBce+3Ec3Z0TGS0jZHmzpU2bw5mPD09UrtLyf7x440Hno3MaVI/mAEAEo0a3igEUfsY9da1fr2mIOubw6w3rVTfXZLJSGeeKX30o9KhQ8EvpMrnpcsvl6r9vre1Seec43zo8DKWRua00bp3AABcUMMbZ41mhquJemMDvzJ3QcxNkMetxC37OdnoqPTGG06WN593PpQEudBr+3Zpxozq9ykWnYV0XrO9jcwpPXEBABEgw5smUfaVJXM3VXn2s5Kw5qdWxrlcUFcH6OwAAPARfXhbUakWNgo9PdLg4PSSilbN3JX3fn3zTad0YWxs4j5hzo/b+yNNtH8rF1T9Nz1xAQARIMMbR0nspxu3zF0zcxjUBg1Rzo/b85dqiHftko4dm1rfG2b9NwAADfKa4SXgjZuoA6NmRFlSUT6ORucwyPmPen4qPX89rzmJH8YAAKlFwJtUUXdbSINm5rBV599LMJ7kD2MAgFSiS0NSxaGfbtI1M4etOv+TeyFX6hjhtRMHO6lVxtwAQCQIeOOGtk3Na2YOkzT/YQdPXj4MsJNaZcwNAESGgDduou6nGzeNBHXNzGFS5j+K4MnLhwF2UquMuQGAyBDwxk3YGyTEWaNBXTNzGOX81xPcRxE8efkw0KolIV4wNwAQGfrwxlGU/XTjpFpQV5qfSl0DmpnD0mNLx16zJviOBOULwgYHpR07KgfbUQRPXnro0o+5MuYGACJDl4akiLodVBTPX2v3tqBbiIXZkaDe7hBu9zdG+ou/kLZt8398XtHJoTLmBgB852uXBmPM88aYIWPMoDGGSDZsUS92KX/+e+6RFi6U1q0LbgyFgvTWW9Nvn5wRC/KyvtuxX3tN2rSpueNWKluolLH94Q/dSxxyOWfjiMmslX7842gXQXkpCWnVTgWUKwFAdKy1Nf9Jel7SOV7ua63VihUrLHx0++3WZjLWOiGN8y+TcW6P6vkla9varJ03z9oXXvD3+V54wTlue/vU52tvn/p8l102fUyStd3dzY+h0rFnzGj89ZZeV2kuM5mJ11NtjsvvW3LzzRM/j+K8aES1OcCE0jlx2WXOV+YHAFxJ6rceYlMWrSVB1Itd3J5fkorFYBZKlbKrY2NTb1+0aGpGLMgWYj09TolAOWsbf73VMtLlC8JKz10sTty3PMN88ODEz0vivgiKTgW1RX1FBwBSyGvAayU9aowZMMbcFuSA4CLq3rBuz18SRIBVKcCePXvq5d8gW4jlclKby69Hsdj46632waX8cvf8+dMfb620fftE4BP1edGIWh/eWrXcYTI+FACA77wGvP+btfZSSaslfc4Yc2X5HYwxtxlj+o0x/ceOHfN1kC0v6t6wped3y3gGEWB5DeSCrInMZqW1a6cHvc283kqva9EiJ7hbs8a5bedO6ROfqJ1hjvq8aES195bMpiPqKzoAkEZe6h4m/5P0f0v6YrX7UMMbgFJNX3d3NDV9L7zg1IzOmFG9rtTrsarVJ8alztPvcbgdb+5c51/5c/zyl85c16pRjvq8qFe9dcyTa5Jbpa416pp9AEgQeazhrdmWzBhzlqQ2a+3I6f/+F0l/Y639p0qPoS1ZipXak1Xqw+rl8V5aMzX7PH7xexzlxxsZkb73vakZvbY26ZxzpHnzpOeec0KekmqtyuodQ9Qt7srntFobup07W6elF+3LAMAzr23JvAS8F0p66PS37ZK+Z639f6s9hoDXR1EHJ36rt99s2lUK8iSpvV06dUqaMcNZwOdH4BPnYKrauSG11nkTlw98ABBzXgPemjutWWsPS+r0ZVSoT727byVBkusTg/jw0dMj/epX0ztSSM5t7e1Oje/s2f4EPl52r4tKLuec3+XBeC7n1Dcn9bxpBLstAoCvaEsWB5VWpqdxtbbfnQXCWtUf1IKqXM7J4FYyNuYEu/m8EwA1G2DH+QNHtUWISexIAQCIjZoZXgSsWha33uAkCeUP1bJ49QojA16a0x/+0OmDWyoB8iMzWjq2W/uzEr+Dup4eZ57KSwPiEji6ZTYLBafWuVh0OldYm4yOFACA2CDDG7VqWdx6slphtXRqNqPqZyuxoDPgk+f06NGpi8dKz9doZnTysd980/0+bW3+B3VJa2VWmqfvfc+pZzbGyYh/6lPJLu0BAISKgDdq1bK49QQnYZQ/+BVUl7J4zV6mrzR3v/iFP2UO5XNarpnMaLVjl4K6T3/a/6AuyN7FlTTyIan0mK4uaXh4Yp6KReeDQEcHwS4AwDNKGqJW7RJzKTjxslo7jNrMuC14qjR3zzzjbLvbbJlDpR3fSs/TTGa00rHPOktavz7YcpQwF0Q1UnZS/phycak5BgAkBhneqNXK4nrNhoaxqCduC57c5q6tzbn07Uem221O29qkc89tPjNa6f1av765rLdbNjXK7Xr7+pz628nvx8hI9fcjyMw6AKAl1ezD2wj68NbJj56bYfRXjWMP3fK5e/xxaf/+6ffr7nYC9nqPHdScBnFst2OeeabzszfeiKbvbmen+/vR2elke91U600cp77BgAdJWEsMJJnXPrxkeOPAj5rWMGoz47jgqXzurrzSv0x3M3NaK6saxPvlVnLy+uvTM6xhtrar9IG6WKz8GLfstzH+ZNaBEIW1lhhAbWR4UZ+47QBVnj5Zu1ZavTrancSi2s2sWmbU7b71ZrwbsXy587rLVcvwxnk3OKAOcbwoBqSNbzutAVPEaQeoSguidu2Stm+PLigPYnGfl+uibov4jHH+Tc6ohlkDe8UV0tNPT91JLpNxbq+knsWaQIzFbdkD0MoIeKNAUZc/KgWW27dHG5T7/VfOa6cDt009KtXwhlWG0uhGI3H6YFUvfr9xWtz3eQFaCTW8Xvm10p2iLv/ENX3id8cMrz2W3eqCh4acf2H23a01pjSXJvD7jUniuOwBaFXU8HrhZ00hRV3+ietc+l2DWqk2N6w6XHgX13MSkYnbsgcgbejS4Cc/dzGLa1YyieKaPvE7qxlGj2X4g99vlPFrY0kAzSHg9cLPP2KtErzk89IllzgB6CWXBJOJjPPl8tJfuZ07ne/XrGm8FCaugX29otwAIyyt8vsNAAlDSYMXfl6mbIWWS/m8dPnlU3uwGiM9+aQTELQKP9/rpF8XbYXzXmqd1wkAMUFJg5/8zLDFOSvpl1tvnb7hgLXO7a3Ez1KYpF8X9XMu4qwVfr8BIIFoS+aF331Bk9xyyYvDh+u7PQhxaA1FPeeEVpqLtP9+A0ACEfB6xR8xh5dA8sILnVZY5S68MLwxeulbGzSacE5gLgAAEaKkAd557TF6//1Oze5kxji3h8Gvy+fNLrJqphQmiQu8qo05LQvvAACJxKI1eFfP4r183qnZPXzYyezef394C9aq9a3dudNbqUM+L111lXTypPN9e7vU0VF/lriRxWZJXPjkZcxJX3gHAIgdr4vWKGmAd/XUYfb0SPv3O/9dCnR6e8Opp610+XzRIm+lDoXC1GBXksbGpJER53XUU9rSSClMtQx1XMtqvIyZsiAAQEQoaYB3jfQYbXar1UYu7Ve6fC55K3Xo65sa7JaMjYWzyCqJC7ySOGYAQMsg4IV3jdRhNlNP22iwXKk11MGD3oKyaptkhLHIKombFyRxzACAlkHAC+8a6THaTOavmWDZrW+t16Csp8ep2S03c2Y4i6ySuMAriWMGALQMAl7Up94NEJrJ/Pl9mdxrUJbLOQvUJge9M2dKu3eHs8gqiZsXJHHMAICWQZcGBKuZjgN+buk8eTxeOgXQUQAAgNjz2qWBgBfBazR4TGJ7LgAAEBrakiE+Gm1H5feWzgAAoCUR8CLe6N0KAACaxKI1AAAApBoBL9CsRjbHAAAAoaGkAWhG+cK6StsVAwCAyJDhBZrhtjnGa69JXV1kewEAiAkCXqAZbptjWCsdPep9K2QAABAoAl6gGW47yZXUsxUyAAAIDAEv0Izy7YrLNbMVMgAA8AUBL9CM0uYYGzZI554rtZX9SmUyTvcGAAAQGQJeoFmlzTH6+6W5cyeyvaWtkHO5aMcHAECLI+AF/DI529vd7XylPRkAAJGjDy/gJ7ZCBgAgdsjwAgAAINUIeAEAAJBqBLwAAABINQJeAAAApBoBLwAAAFKNgBcAAACpRsALAACAVCPgBQAAQKoR8AIAACDVCHgBAACQagS8AAAASDUCXgAAAKQaAS8AIBUKBam3V+rudr4WClGPCEBctEc9AAAAmlUoSJ2d0okT0uioNDgo7dgh7dsnZbNRjw5A1MjwAgASr69vItiVnK8nTji3AwABLwAg8fL5iWC3ZHRU2rMnmvEAiBcCXgBA4vX0SJnM1NsyGaeeFwAIeAEAiZfLSXPmTAS9mYzzfS4X7bgAxAMBL1oPS7mB1MlmnQVqGzY4v9obNrBgDcAEY631/aBdXV22v7/f9+MCTStfyl1KA/GXEQCAxDHGDFhru2rdjwwvWgtLuQEAaDmeA15jzAxjzK+MMf8Y5ICAQLGUGwCAllNPhvfzkg4FNRAgFCzlBgCg5XgKeI0xF0i6VtLfBzscIGAs5QYAoOV4zfDeLSknqVjpDsaY24wx/caY/mPHjvkyOMB3LOUGAKDltNe6gzHmOklHrbUDxphVle5nrb1P0n2S06XBtxECfstmpS1boh4FAAAIiZcM70pJ1xtjnpf0fUlXG2O2BzoqAAAAwCc1A15r7X+11l5grV0g6UZJP7PWrg18ZAAAAIAP6MMLAACAVKtZwzuZtfbnkn4eyEgAAACAAJDhBQAAQKoR8AIAACDVCHgBAACQagS8AAAASDUCXgAAAKQaAS8AAABSjYAXAAAAqUbACwAAgFQj4AUAAECqEfACAAAg1Qh4AQAAkGoEvAAAAEg1Al4AAACkGgEvAAAAUo2AFwAAAKlGwAsAAIBUI+AFAABAqhHwAgAAINUIeAEAAJBqBLwAAABINQJeAAAApBoBLwAAAFKNgBcAAACpRsALAACAVCPgBQAAQKoR8AIAACDVCHgBAACQagS8AAAASDUCXgAAAJ8UClJvr9Td7XwtFKIeESSpPeoBAAAApEGhIHV2SidOSKOj0sCA9K1vSWvXSps3S9ls1CNsXWTPc0LTAAAgAElEQVR4AQAAfNDXNxHsSlKxKJ06JT3wgBMIk+2NDgEvAACAD/L5iWB3MmudQLivL/wxwUHACwAA4IOeHimTcf/Z6Ki0Z0+448EEAl4AAIDTmll0lstJc+ZIbS7RVSbjHBPRIOAFAADQxKKzrVulvXudr/XU3maz0r590qc/Lc2YIRnj3J7JOIFwLhfc2FEdAS8AAAhcEtp1lS86Gx11r72t9lqyWWnbNum3v5U+9znnPhs2OIEwXRqiY6y1vh+0q6vL9vf3+35cAACQPOXtukoZz7gFgd3dTmbX7fZ83vnvpLyWVmGMGbDWdtW6HxleAABiLO6ZUS/j85o5jZrborPy2tukvBZMRYYXAIAGFApOkJPPO4HS2rXS9u0T3+dyzWf84p5N9Do+L5nTOPDyepLyWloFGV4AQCDinnEMQ/nipnvvlS6/3PnayGKnSuKeTfQ6Pi+Z0zgoLTrbsKFy7W1SXgumIsMLAPAs7hnHsPT2OkGt2yYDJZmMEzBt2dL488Qhm1ieyZ6cufY6vjSdN2l6LWlAhhcA4Lu4ZxzDUmlHrcm8bjRQLWMeRDaxngx9KbgrZa6/+U3pT/90Ipj1Oj4vmdOkSNNraSVkeAEAnsUh4xgHfmV4a2UL/c4m1nu83l4n2B0bm3r7zJnS//pfzn+T7USUyPACAHxH/aKjtKNWaS4yGWeTgfb2ie+9bDRQK2Pudzax3gx9Pj892JWkkyedx5DtRFK0Rz0AAEBy5HLSjh3TM3qttoNUKdDr63PKFhYtkv7wB+nxx52fr14tbd5cO/BzK42YXApRrX62EbWer1xPj3tGX5p4TDbbXJ0yEAYyvAAAz8joTSgFejt3Sj/5ifTww9LRo9Lvf+9870W1jHmz29zW+3xucjmnfKFce3vrZfWRbAS8AIC6lAK9fN752orB7mTNLORzK40oZcyDWCBY7fncZLPS7t1Tg972dqmjo/Wy+kg2Al4AAJpQb5nAZNUy5s0ct5Hnq6Snx1mgdvvtzmM2bgwvq0/PZ/iFGl4AANR4vWxPjzQ4ODU4rWchX6Ua2MWL3etnFy3ydtx6n8/vxzSrvKPE4KBTP14p2Pa73hnpQsALAGh59QZXk7GQLxhuJR2//7107bXST3869X1p5v1Da6CkAQDQ8pqplw1qId/Bg+63HzrU3HFrqbeMIKiyg0qbewwNTV+8x4YoqIUMLwCg5TVbLxvEJf9mSyUa0UgZQSOZVS/lB26vv6QUzJbmPIh657ijhKM+ZHgBAC1v8WL325utl62lWna03o4Kfqg3U9pIZtVru7XS63dTHsy22oYoQbSsSzsCXgAA6uDXJfxaQUsUPY/rzZQ2kln1GiSXXv+yZdOPUR7MRvHhIEqUcNSPkgYAQMvzWi/r1+KoQsFZfPX730/cNjloKV2qD7s7Qr1lFI2UXdQTJGezzgK1yXPuFsyW73zX3Z3uS/ytWMLRLDK8AICW5/WSuB+ZtVLQPDQ0/WdRBy31ZkobyazWW37gNdPdShuitFoJhx+Mtdb3g3Z1ddn+/n7fjwsAQBDKM7elwK08sOrudu+N293tBFpe9PY65Qtui7GMkf7iL6Rt2xp6Gb4oLYbymilt5P5e5hqVMYcTjDED1tquWvejpAEAkEr1rGL3ekncj84JldptSZK10o9/7Iw9qsCl3jKKeu5fek+yWee1trVJV1yR7vKDILRaCYcfyPACQAtplVZGQWXAvB632jxXy/BKzjE3bAh/Z7Oglc9de7s0Y4Z00UXSlVdWPhdb5ZxFY7xmeAl4AaBFtNJlULeg0i2QLBSkTZukXbuc71evljZvbu4Sfq15zuelq66STp6s/Bz1lEgkRbVAv9oHh1Y5Z9EYrwEvi9YAoEW0UisjL6vYCwWn5dV3vysdPer8++53nduqtRqrtTiq2jwXCk5QfepU5eN7LZEIaoezoFQr5ah0LrbSOYtgEfACQItopVZGXlax9/VJIyPTH/v6680FVNXmuRTAjY1N/Xlb28QYvfSPTeLGA27vyWSlOZocyP/wh61zziJYNQNeY8wsY8weY8w+Y8zTxpj/J4yBAQD81UqtjLy0y8rnpWJx+mOtbS6gqjbPlbKc55wzveVWtQxuEjOfpfekvcJy+fZ2Z2e7yYH88ePT75fWcxbB8pLhPSnpamttp6Tlkj5ijHlfsMMCAPgt6btR1XMJ30vv1p6eiczqZMY0F1DlctKZZ04cu63N+T6XqxwMf+ITU0skamVwk5itL70nlbZrnjHD+To5kC99IDHG+VrvOZu0sg8Ep2bAax0nTn+bOf3P/5VuAIBARbFVrV8auYRfq9Y2l5M6OqY/7h3v8OdDQGlNuLXSm29K//E/OiUUZ55Z+0NHrQxuUrP12aw0a5b7zy6+2Nnxzi0DPn9+/eesn2UfBM7J56mG1xgzwxgzKOmopH+x1k5bO2qMuc0Y02+M6T927Jjf4wQA+CCpu1EFcQk/m3V2O7v5Zuncc51/N9/s3NbMvPT1SW+8MTXgffttJ1B74AEn6L3wQif4qhTA1crgVsvWxz04qxSsX3GF9wy4F36dM0msl8Z0dbUlM8a8U9JDknqttQcq3Y+2ZAAAP/mxw1lYKo11srY2ae7cytlKL23V3NqjSfFv41Wt1Zjk3/j9Ome8trhDNAJpS2atfU3SY5I+0ujAAACol1vmr62tcj1olGp1I5Cc2tRq2UYv9dZu2fokLGarVlrjZ9mNX2UfSayXxnQ1M7zGmPmSRq21rxljZkt6VNJd1tp/rPQYMrwAAD+VeuYOD0+9fe7c5ksQ/FaewaymWrax1gYXlY6XlEx40PzatIIMb7z5meE9T9Jjxpj9kvbKqeGtGOwCAMIR91rNZpS/Nkn66EcnVuuXvPFGvLKX0tQsZWenNHPm9HFLtbONpQzuzp3O92vW1H6fw1jMlpTzzq9scdK7m8DB1sIAkEBp3nK10mvLZqX9+6ffP+7Zy9L2xdu3OwvYikXnNbW1SRddJF15ZeXsbb3vc9DnRZDHL2W083kncPeS0W7kMc2MrZ5sO8LB1sIAkGJJqNVsVKXXZm1yW3Ft2yb99rfSZz/rBIxtbc72wvv3V1/1X+/7HHTruaDOu0Y6IYTZPSGq7iZJyaYnAQEvACRQmhfSVHptbW3JvrRcCpquuMLJ8pa2F64WNFaai1/8ovbzlIIzyb+gKajzrpFAOs0f+iTaofmNgBcAEiipGw94Ua1Pa1I3zpisnqCxp8d9K95nnvEW+PgdNAV13jUSSCfxQ189Gdu0B/RhI+AFgARK80Kaaq8tqRtnTFZP0JjLTWy5O9mpU94CH7+DpqDOu0YC6aR96Kv3w0cSA/o4I+AFgAQq1Wp+6lPODmHz5knXXx/1qPwR5hbI5Rm3fD74msl6gsZs1lnYVm5szD3wKX89v/hF9aCp3hrRoN6bRgLppH3oq/fDR9IC+rijSwMAJFRUnRrCWhkftPL5a293MqczZjgBZZDzWc+qf699YN3Oh9LiuFK98OTH5nLx6vTRSCeEJHVPqLdHcpo7sfjJa5cGWWt9/7dixQoLAAjW7bdbm8lY6/QvcP5lMs7tQXnhBWvnzZt43kzG+f6FFyZ+fvvt1l52mfO1dPvkx1f7eZjc5q/8X9Dz6UWtOS9xez3t7dbOnOn+2Eqvf9myaN+XtGrk97X0PnV3R//7EleS+q2H2JQMLwAkVBS7alXLNtbKGMYtY1Vp/tzuF3WfXy+ZzEqvp7PTWfBX/thqr3/uXGejj4MHk53Fj5O4nf9pQR9eAEi5KGr8qi2kqVWj6NcCKr96k7rNX7m41ExWWqw3eS7eemt6R4dSdwu3x1Z7/cPD0oMP0g7LTaPnX5i16ZiODC8AJFQUGaNqGd58vnrG2Y+MtJ+v2e1YY2Ph1PD6odka5NLjf//72s/lVjPc7NiTWAdOljZ+yPACQMpFkTGqtjK+VsbZj4y0n2223ObvySeljRuTkYErn4uxMacy9D3v8Tb+0ut36wJRrp52WLUyoEnaUKH8tWzaRG/cpCLDCwARKAwX1PdEn/JH8uo5v0e5lTll58Y0sipTqZ60VvbLj+xYFHXLcVWtBveXv3Q+YNRSKEjLljklDNV4zfB6eY+9dp2ImttrKRadLHq5Vjz/4oIMLwDEVGG4oM57O7V1YKv2vrRXWwe2qvPeThWGY5jiclGpnrRWxtmPjHSSepP6VWtcSbWA9tZbvR2jr096442ptxkjnXFGY/1tvWTgk7KhgttrsdaZn8niev5hKjK8AFpCnDKqvY/0auvAVo0WJ/7qZ9oy2rBig7ZcE6MUVwwlpYYyjHEWCk75gpv2dunw4drPVW9nh1q8ZOCTkuGt9FpmzHD6G8f5/GslZHgB4LS4ZVTzR/JTgl1JGi2Oas+RmKW4YigpK9393tLXTaVd2CSnnrdWXWyh4HR2KFets8Pkx7plr71k4JOyQ1ql17J2bfzPP0xHhhdA6sUtoxq38cB/YdUa5/PS5Zc7l9rLVcualjLQIyPTd2GrlbGslr2WvGW2k7BDWlKuJrQ6MrwAcFrcMqq5lTnNOWOOMm1O+ijTltGcM+YotzJmKS40LKxa454ep7PErFnTf1atLraUgZ4c7ErSxRfXDnavvdZpZeaWvfaaga9UBx4nSbmaAG/aa98FAJKt5/weDb48OC2j2n1+NCtNsnOz2rdxn/qe6NOeI3vUfX53oro0tDovPWRzOWnHjunZwSAu2/f0SLfc4l4XWynAdls4JkmzZ9fO7Lr17Z0cXJeC2TRI02tpdQS8AFIvtzKnHUM7dOLtExotjsYio5qdm6V8IYHKL3MPDjqBbXnmr5QdDOuyfaUAe+1ap8a2PDjv6XHG7jVAliaywm7oVIC4o4YXQEsodWkgo4pmeOkwENUuYuV1sWvXSqtXN1drO1m1vr/z5nG5H9HwWsNLhhdASyCjCj/U6iHrNQMchPLL7729lTtFbNlSfwbaLSssORtX/PSnBLuINxatAQDgUa3FaGG0I/OqVnBe78Ixt3Zi8+YR7CIZCHgBAPCoVg/ZOO0i5neniGxW2rXL6eRw1lnO1127CHaRDAS8AAB4VKtVVZy2PvZ7g4dCwakJfuYZ6Q9/cL6uXj2x6UTQWykDzWDRGgAAPonbZgV+bvBQbcFeLhev143WwaI1AABCFnY7Mi/j8auPbLVyjWq1y/SxRRwQ8AIA4KO0blZQrXdvnGqXATfU8AIAgJqq1QTHqXYZcEPACwAAaqq2YM/vBXKA3yhpAAAAnlQq14hb7TJQjoAXAIDTotoWOA3SWruMdCDgBQBA0W4LDCBY1PACAKB4bQvcCtioItmS9v6R4QUAQLTWChPZ9GRL4vtHhhcAANFaK0xk05Mtie8fAS8AAKK1VpjIpidbEt8/Al4AAFS9zyz8RTY92ZL4/hlrre8H7erqsv39/b4fFwAAJF95DWgpm84HjGSI0/tnjBmw1nbVuh8ZXgAAECqy6cmWxPePDC8AAAASiQwvAAAAIAJeAAAApBwBLwAAAFKNgBcAgEmStmUqgNrYWhgAgNOSuGWqnwoFZ7esfN7ptZrLtcbrRvqR4QUA4LQkbpnql1Kwv3WrtHev87Wzkww30oGAFwCA05K4ZapfWjnYR/oR8AIAEs+vutskbpnql1YO9pF+BLwAgETz81J8LudskVoKektbpuZy/o45jlo52Ef6EfACABLNz0vxSdwy1S+tHOwj/ejSAABINL8vxWez0pYtzY8raUrBfl+fM3fd3XRpQHoQ8AIAEq2nx2kfNjno5VJ8Y1o12Ef6UdIAAEg0LsUDqIWAFwCQaK1cdwvAG0oaAACJx6V4ANWQ4QUAAECqEfACAAAg1Qh4AQAAkGoEvAAAAEg1Al4AAACkGgEvAAAAUq1mwGuMyRpjHjPGHDTGPG2M+XwYAwMAAAD84KUP75ikL1hrnzLGdEgaMMb8i7X2YMBjAwAAAJpWM8Nrrf2dtfap0/89IumQpPODHhgAAADgh7pqeI0xCyT9uaS8y89uM8b0G2P6jx075s/oAAAAgCZ5DniNMXMk/UjSndba18t/bq29z1rbZa3tmj9/vp9jBAAAABrmKeA1xmTkBLs7rLX/EOyQAAAAAP946dJgJH1b0iFr7d8GPyQAAADAP14yvCslfVrS1caYwdP/rgl4XAAAAIAvarYls9b+myQTwlgAAAAA37HTGgAAAFKNgBcAAACpRsALAACAVCPgBQAAQKoR8AIAACDVCHgBAACQagS8AAAASDUCXgAAAKQaAS8AAABSjYAXQOIUhgvqfaRX3fd3q/eRXhWGC1EPCQAQYzW3FgaAOCkMF9R5b6dOvH1Co8VRDb48qB1DO7Rv4z5l52ajHh4AIIbI8AJIlL4n+saDXUkaLY7qxNsn1PdEX8QjAwDEFQEvgETJH8mPB7slo8VR7TmyJ6IRAQDijoAXQKL0nN+jTFtmym2Ztoy6z++OaEQAgLgj4AWQKLmVOc05Y8540Jtpy2jOGXOUW5mLeGQAgLgi4AWQKNm5We3buE8bVmxQ97u7tWHFBhasAQCqoksDgMTJzs1qyzVboh4GACAhyPACAAAg1Qh4AaAKNrkAgOSjpAEAKmCTCwBIBzK8AFABm1wAQDoQ8AJABWxyAQDpQMALIBJJqI1lkwsASAdjrfX9oF1dXba/v9/34wJIh/La2NLmEXGrjU3KOAGgVRljBqy1XbXuR4YXQOiSUhvLJhcAkA50aQAQuiTVxrLJBQAkHxleAKGjNhYAECYCXgChy63Mac4Zc8aD3lJtbG5lLuKRAQDSiIAXQOiojQUAhIkaXgCRaOXa2MJwQX1P9Cl/JK+e83uUW5kj2AeAABHwAoAPvAaxbFcMAOGjpAEAmlQKYrcObNXel/Zq68BWdd7b6bqZRlJasgFAmhDwAgkT9g5lSdgRLWr1BLFJaskGAGlBSQOQIGFfDo/68ntSal3rCWJ7zu/R4MuDU+5PSzYACBYZXiBBwr4cHuXl93rKBKJWT19hWrIBQPgIeIEECftyeJSX3ysF25t+til2JRb1BLG0ZAOA8FHSACRI2JfDo7z8XinY3j60XW2mLVYdDkpBbN8TfdpzZI+6z++uWn7Ryi3ZACAKZHiRKmlfYBX25XC/n6+e98etTMDIyFobyw4HpSA2f2teW67ZQsYWAGLEWGt9P2hXV5ft7+/3/bhANeULrErBWdTZP7+VFnJ5ySTG6fnqfX/c7l+0RZ2yp6bdt/vd3crfmm/o9QEAkssYM2Ct7ap1P0oakBrVFlgFcfk4qg4CYV8O9+v5ar0/bvNZXiYw8vaIvjf0PTocAADqQsCL1AhzgZUf7bqS0nLLL9Xen2rzOTnYLgwX9JNnfzItS0yHAwBANdTwIjXqaQ3VrGbbdSWp5ZZfqr0/bvM5/Nawuu7rmlLrS4cDAEAjCHiRGmEu6Go2m9yK28tWe3/c5rOooo6+cXTah4HJi8NyK3Pqe6JPy+9drku+dYk67+30bbGi3wsg076gEgDijJIGpEa9raGa0Wy7rlbcXrba+7N4/mLtfWmv6+Mq1WKXsuQjb49orDg2fvvBowebblXm9w5zUe9YBwCtjgwvUiWs1lDNZpPDLL+Ik4rvT41mMW4fBkpZ8snBriSN2bGms+V+Z+BbMaMPAHFCwAs0oNlaUj/KL9J0ifzg8YNVf+72YcAtS17SbLbc7wy8X8dL03sOAGGipAFoUDPtupotv0jbJXK3EpGSSh8Gaj2mmWx5pZKVN8feVPf93XV31fBjx7q0vecAECY2ngCaEFVrsd5HerV1YOu0AGrDig3TgvAktD9z22SizbTp4nMu1hXvucJ1zJVqeNtNuzpmdvhaw5tpy2isOKYZZobG7Fjdm5r4sSlKPe85ALQKrxtPEPACDYpyZ7fu+7tdF3mV7zgW993nJgfji89ZLBnp0LFDnjPepcf/4oVfqGiLajNtFQPkRse258gevTn2pg4dO6QxOxFY1xtsNrtjndf3HABaCTutAZMEkeUMe2e3ybxeIo9yjLW4XaKvNxgPcte5ycfuvr97SrAr1V+D2+xY6ymLSEJWHwDCxKI1pF5Qmzx4WYgU1CIjr4ve4tz+LEmdC/zoqtHsueD1PW/FTU0AoBYCXlSVhlXhQQVWtYKgIAMPr10i4tz+LIxg3K/zt9muGl7OhVpj9fqeJ+mDBACEhRpeVBT3+k+vgqp9rDU/bouM2tSmjy36mP7DnP8QyuXmOL+HjSzCqudSvd+vvZka3Fqv1c+xUusLoJVQw4umxbn+sx5+tIRyU6u1WKXtcnce2ql2064xOxZ4a6kwd5+rV25lTjuGdkwL8iplTetty1Xp/N302CZ1nNHhOWj2oxa2Vjbbz9+1oM53AEgyAl5UFOf6z1qmrP6fv1hnZs7UG6NveAqs6lFtIVLP+T0aeGlARRWn/ay0AKrZDxFeArIgF3Y1o95gvN6gsNL5u33/drWZtmlBc+k5SnO59pK1Wr1jtS99b2sFoX7+rpV/kDAyKtqiRt4eUWG4EIsPOwAQNgJeVJSUTFF50OcWqJyZOVOfWvopHTruveVVs3Irc/pW/7ca2jbXizA2Igh6tX89wXi9QaHb+dumNhVtUafsqfHHv/bWa/rrf/5r/c/f/s8pc3n/U/frlD013uO3mQ8ntbLZfv6ulT5IbPrZJm0f2i5rrU7ZU/re0Pf0k2d/EotyFgAIG4vWUJEf2982q9ZCHrfFQFdtu0ojb49MyQS+MfqGOmZ2KH9rXluu2RLKH/zs3KzWLlsrI1P1fo0GNkEvTgprtb/XhWWL5y9WW9n/sqrNndv5a4yRLfsEYmX1o0M/0sjJqefMyVMnp2xoUbq9kQ8ntRac+f27lp2bVcfMDrWZtvErDCxeA9DKyPCioqjrP71kMN2CPjdRlWJsvnqzfvLcT2ru2NVIYBN0yUkYNdxes9SF4YJ+/MyPp5WHnJk5s+LcuZ2/IydH9MD+B1yD3vI+u26aucJRLZsdxO9akkuSAMBvBLyoKsr6Ty8Bl9sfdTdRlWK4BTJrL1mr7fu3Nx3YBF1yEkbA5DWo7nuiT2+MvjHlsUZGH73oo1Xnrvz8LQwXtH1o+3hJQzXtpl0z2maoaIu+1357GWuzklKSBABhIOBFbHkJuNz+qIcdqNTiFsj0XNDT9HHr7XJQrzACJq9Btdv9rKwOHT9U1/Nl52a19pK1enDfg1OyxZm2jHP5v+yc2XXTLl8+nEQh6PMDAJKEGl7ElpdNE9xqHztmdmj3ut01G/R7EeeNN7xuRNCoZupKvc6bl/e4MFzQW2NvTXts6X71vkeb//fNmjtr7rTX5XbO9FzQoy3XbAm19tsvlc4PSbE9pwEgKGw8gdjy2oy/mQ0B/Hj+NGtkbkvzNvL2yPiiLyOjNRev0d9+5G+nvXfV5tjtWJKTxe+Y2aFdN+2a0pHD63sU1DlTS9BdL7w8f6uf0wDSxevGEzUDXmPMdyRdJ+motXaplycn4IVfogpMpMZ2Aks7LwFb7yO9unfg3mkdDiRp7sy5GvqrIc8fWNzeA0ladu4y/fRTP1XfE32JeY/iEGxyTgNIGz93Wtsm6ZuSHmh2UEC9olw0l5ZV7m5BqqS6M41eOyrkj+Rdg11Jev3k69MWpLktLOt9pFf5I3n9+/C/uy5KnN0+W9m52US9R3HYuTBJ8wUAfqoZ8FprHzfGLAh+KEC8pGGVu1uQ+uD+ByVpfOc5rxtWeA3Yes7v0d6X9roew8pWDa7Kx+vWw3jye5Ck9ygOwWaS5gsA/OTbojVjzG3GmH5jTP+xY8f8OiwQmSA23gh7EZxbkDp8clivn3y97g0rvAZsuZU5zZwx0/UYRqZqcFU+3lK/3NKGE+XvQRw2R/HKywK9oCVpvgDAT74FvNba+6y1Xdbarvnz5/t1WCAyfndBCGvnsskq9Sku33jBS6bRa8CWnZvV7nW7dUbbGdOO8Y6Z76gaXFUa7zlnnuP6HgTdqcJPcQg2kzRfAOAnT10aTpc0/COL1oDGRbFgqPeRXv23vf9t2g5l5byMo95FV4XhgjY9tkm7fr1LstLqP1mtzVdvrhpcpX1RVZSLMAEgjXzr0nD6YAtEwAs0pfv+btfa1u53dyt/az6Q5ywMF7TwGwur7ixWT7eAoAO2OHQyAAAkh29dGowx/0PSKknnGGNelPR/WWu/3fwQgdYSxYKh7Nys1i5bqwf2PzCljKHdtGvR/EWa3T67rsC1nq4ZjfScdduKmSyod1H3+QWAuGLjCSAkUWUvo3heMrXhY84BtCKvGV62FgZCEtWCIbfn3XXTLvU90RdYt4hqLczgXT1dPZhzAKiMDC+QcPVexg4jExhFvXKj4loGUO/7lKQ5BwC/kOEFWkAjrc7CyATGoeesF1G0ivOq3vcpKXMOAFEg4AUSrJHgNYwdv+LQc9aLqMoACsMFrfvxOr3ra+/Su/6/d2ndQ+umBdn1vk9JmXMAiAIBL5BgjQSvYWQCk7LBQRTb/RaGC1r2rWX67r7v6ugfjuroG0f13f3f1bJvLZsS9Nb7PiVlzgEgCjXbkgHwn191o420OsutzGnH0I5ptaHNbpns9nrivllEFK3i+p7o08jJkWm3v37ydfU90Tc+Z428T0mYcwCIAovWgJCUgsLHX3hczx5/VqfsKY0Vx5paNFa+sKlNbTLGaO2ytVV3NfNzA4kkt8MKeuxuHwTW/HCN6+IyafoCM3ZmA4DqfN1prV4EvK0rriveo1YeWJVrZvvc0ha+2/dvV9EWZWVDDTqTvh1wUEFlpWD6+ouu14P7Hpy23bOR0ecu+1wi5gwA4sK3ndYAr8r/wA++PKgdQzsSkekLWvniqHLN1I1m52bVcUaH2kzb+BbCkxdfBR1ARVEH66egygAqLYiTlTpmdmj45PCU+79j5jtYYAYAAWHRGnxD4/vK3ILCyZqtG40y6Fw8f7Hr7YvmLwr8ueOs0nty6FCYISMAAAlISURBVPghDf3VkG7uvFnnnnWuzj3zXN18yc0a+quhlv9gCABBIcML3yQ90xckt8VRJX4sGoti8dW4SlVR/ldL+Sro8ptq70l2blbbPrrNt+cCAFRHhhe+ofF9ZeU9UttNu2bOmKnOd3X60j4qyh6sB48fdL390PFDvhy/nu116zlm0BtO0BcXAOKDDC98E0S7q7Qo9UgNasV90MevJsjsclB14dXKb/yq543yPQEATEWXBvgqiW2U6CzRHL9ae7m9D31P9AXSAaL7/m7X1mDlbcEAAPFGlwZEImmN7+ks0Tw/Mpml3cdeP/m6rKz6X+rXg/sf1B+/848DqQuPtOYZABA6Al60tDAubbeCZj/obPrZpiltuqyshk8O6+TYSWXaMr4HppTfAEBrYdEaWlqcO0sEsVgrrnb9epfr7cffOB7Iwq9SVnrDig3qfne3LwsHAQDxRYYXLS2ul7ZbrtTCuN88o22GfnXrrwKpC09a+Q0AoHFkeNHS4to6Ko2beFTLWK/+k9Wuj1n9J6vHA9P8rXltuWZLOgP+OsU9+x/38QFoPXRpQMuLY2eJtHURqNXJobRobeTkiIoqqk1t6pjZwe5jLvzqitGq4wOQLl67NJDhRcuLYwYxbZt41MpYZ+dmNfRXQ/rsZZ9V97u79dnLPkuwW0Hcs/9xHx+A1kQNLxBDaesi4GVxIDW13sR5oaUU//EBaE1keIEYSlsXgbRlrKMU97mM+/gAtCZqeAEEjrpO/8R9LuM+PgDpQg0vgNhIW8Y6aNW6HPg9l353VOC9BhBHZHgBIEbCzJCSjQWQdGR4ASCBwuxyQEcFAK2CgBfANGwcEJ0wuxzQUQFAq6AtGYApWm5b45gJc7vruG6tDQB+I8MLYAouc9cWZAY8zO2u47q1NgD4jQwv4KPSNsX5I3n1nN8Ti22K68Vl7uqCzoCXuhzU2u7aj3PN63MBQNLRpQHwSVpWvPc+0qutA1unXebesGJDondC8+vDSBzmJy3nGgA0iy4NQMjSUgqQxsvcpQBx68BW7X1pr7YObFXnvZ0NlSLEIQOelnMNAMJCwAv4pJ5AKM5dENK4cYCfAWIcts6NQ9ANAElCDS/gE68r3pPQBSE7N5vo8oVyfgaIuZU57RjaMa2cIMwMON0VAKA+ZHgBn3gtBeBydPj8zMrGIQOexrITAAgSi9YAH5UWRlVb8d59f7f2vrR32mO7392t/K35sIbaUtK4yMvLuYbK0tBRBYD3RWsEvEDI4rDKv5I0BwEEiChJ4wcgoFUR8AIxFdc/tnEdF+C3OH/oBFAf2pIBMRWHGlA3SaktjnOHCyQDXS6A1kOXBiACceyCkIQgIAkdLhB/dLkAWg8ZXgCS4tFftpakZKERb3S5AFoPAS8ASckIApKQhUb8xbWsCEBwKGkAIGkiCPDaySCKjg5cioZf4lhWBCA4dGlARWluUYXmRNXRgU4SAIDJ6NKAppQCi60DW7X3pb3aOrBVnfd2siIekqKrpeVSNACgEZQ0wFW1gIbLgIiylpZL0QCAepHhhSsWB6GaJHR0SDN6EQNAfQh44YqABtUkoaNDWlFuBAD1I+CFKwIaVEMtbXToRQwA9aOGF67qbVGF1kMtbTQoNwKA+hHwoiICGiB+6EUMAPWjpAEAEoRyIwCoHwEvACQI9dMAUD9KGgAgYSg3AoD6kOEFAABAqhHwAgAAINUIeAEAAJBqBLwAAABINQJeAAAApBoBLwAAAFLNU8BrjPmIMeZZY8yvjTH/JehBAQAAAH6pGfAaY2ZIukfSakmLJf1nY8zioAcGAAAA+MFLhrdb0q+ttYettW9L+r6k/xTssAAAAAB/eAl4z5dUmPT9i6dvAwAAAGLPt0VrxpjbjDH9xpj+Y8eO+XVYAAAAoCleAt4jkrKTvr/g9G1TWGvvs9Z2WWu75s+f79f4AAAAgKZ4CXj3SvpTY8xCY8wZkm6U9JNghwUAAAD4o73WHay1Y8aY2yX9s6QZkr5jrX068JEBAAAAPqgZ8EqStfYRSY8EPBYAAADAd+y0BgAAgFQz1lr/D2rMMUn/7vuBw3OOpONRDwKJxfmDRnHuoBmcP2hUks+dP7bW1uyWEEjAm3TGmH5rbVfU40Aycf6gUZw7aAbnDxrVCucOJQ0AAABINQJeAAAApBoBr7v7oh4AEo3zB43i3EEzOH/QqNSfO9TwAgAAINXI8AIAACDVCHgBAACQagS8kxhjPmKMedYY82tjzH+JejxIDmNM1hjzmDHmoDHmaWPM56MeE5LFGDPDGPMrY8w/Rj0WJIsx5p3GmJ3GmGeMMYeMMZdHPSYkhzHm/zz9d+uAMeZ/GGNmRT2mIBDwnmaMmSHpHkmrJS2W9J+NMYujHRUSZEzSF6y1iyW9T9LnOH9Qp89LOhT1IJBI35D0T9baiyV1ivMIHhljzpd0h6Qua+1SSTMk3RjtqIJBwDuhW9KvrbWHrbVvS/q+pP8U8ZiQENba31lrnzr93yNy/uCcH+2okBTGmAskXSvp76MeC5LFGDNX0pWSvi1J1tq3rbWvRTsqJEy7pNnGmHZJZ0p6KeLxBIKAd8L5kgqTvn9RBCxogDFmgaQ/l5SPdiRIkLsl5SQVox4IEmehpGOS/vvpkpi/N8acFfWgkAzW2iOSvibpBUm/kzRsrX002lEFg4AX8JExZo6kH0m601r7etTjQfwZY66TdNRaOxD1WJBI7ZIulfQta+2fS/qDJNagwBNjzDw5V7MXSnq3pLOMMWujHVUwCHgnHJGUnfT9BadvAzwxxmTkBLs7rLX/EPV4kBgrJV1vjHleTinV1caY7dEOCQnyoqQXrbWlK0o75QTAgBf/h6TfWmuPWWtHJf2DpPdHPKZAEPBO2CvpT40xC40xZ8gp2v5JxGNCQhhjjJwaukPW2r+NejxIDmvtf7XWXmCtXSDn/zs/s9amMsMC/1lrX5ZUMMZcdPqmD0g6GOGQkCwvSHqfMebM03/HPqCULnpsj3oAcWGtHTPG3C7pn+WsUvyOtfbpiIeF5Fgp6dOShowxg6dv+5K19pEIxwSgNfRK2nE6WXNY0vqIx4OEsNbmjTE7JT0lp9vQr5TSbYbZWhgAAACpRkkDAAAAUo2AFwAAAKlGwAsAAIBUI+AFAABAqhHwAgAAINUIeAEAAJBqBLwAAABItf8fSscJyrsLFKsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa78ef54c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualization\n",
    "cluster1 = X[np.where(idx == 0)[0],:]\n",
    "cluster2 = X[np.where(idx == 1)[0],:]\n",
    "cluster3 = X[np.where(idx == 2)[0],:]\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter(cluster1[:,0], cluster1[:,1], s=30, color='r', label='Cluster 1')\n",
    "ax.scatter(cluster2[:,0], cluster2[:,1], s=30, color='g', label='Cluster 2')\n",
    "ax.scatter(cluster3[:,0], cluster3[:,1], s=30, color='b', label='Cluster 3')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bonus:\n",
    "    Enlist 5 different prantical applications of k means.\n",
    "    Implement Group D excersice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Group B\n",
    "Our next task is to apply K-means to image compression. The intuition here is that we can use clustering to find a \n",
    "small number of colors that are most representative of the image, and map the original 24-bit colors to a \n",
    "lower-dimensional color space using the cluster assignments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADaCAMAAABqzqVhAAAB41BMVEX/////zJn//wDM////mcwA/wD/AAAAAAD/nND/0Jz/0p3G+fl6Wzt/enqZUXXs7+7CwsI+Pj719fXT///k5ADa2tpWVlYlJSW5kWhxWkT5+vxxRFqTk5NacXHLy8tZZmaIk5OxZYuYm5nk4+OLi4tjYz1kZEhKShWlpZJ6iooA+AAA8AARERFLS0vmirj/oddeXl6+AAAA5QDPzwCCoqIA1wAAVACkpKSbwsJycgAAtgDu7gDa2gBEVVXWgKueAAAAjgDtAAAAzAC86+sAagDIeKDKAACiYoKvcpCvkHKigmJgYAC1bZGVWXcAUQCr1tZqamo1NTV7ewAAvwBsAACwsrHHqKiLYmIAqgAAYACKinRaNkhSUgBWVioAfACrAABHKzmQc1ZmPVIAQgAAhgDPpXyZmQC1tQBlf398SmNvi4sfHwA9TU0AcwB4eGB3AABHAAAALwCDAABVVQCQkAC9vQBWRDSti4smFx4AmwAnSyc9nD09Zz1rYWtmFBSKLy8AFwBkQ0NLPEtvsG8udy52nHYyAB6Fk4XFzsXOxM5XuVcnMicfBwA6OgC+vnN6enJ9fT8/MiWtaGg5AACSGhprMTEsLAChoU96ejFXalc4cziYTEyVypWksKQIIwhMGhooYygcHpgzAAAU4klEQVR4nOVdi18iV5ZWHlI26/ISVBDoiEk6oKW8HwoqZUlGCvFFAsEo4BMi07atTrtJemfz2Elv0tuT7ZnMpmc2u3/qnlvFU5CHQBXo9/vdoupwkftxzrn3u9d69PU1hMeffPrpJ6j84dM/0OVTpaGxT/YUnHtLN2H+q4TrVrUdZ0sDVbAwyXW72gwn0DQ/fboysFTkuwLl7/fLoxEz0Ho7sPR2ZQH2GKoLC7BZest109qKT4DayhnQWjo7W1rY3l4yn2+f0aH8lOumtROSc9qBZ2fmgQXz0gK4cu984Bz5c8Ac4bpxbQTDc2BpZdu8Zx7Ye/pyAV4YnitSrhvXRtA8zQuoLJjN5wN7JTzv0xiq/i/kzZfne2dL5u297b2X29Admen8XLhP/uyLrNBhu7KEOlt4oUcXuts947ppbYXyZTWVANgb4rpp7YVhryrNpXtGEzz6lyrC7755E8HwdnuBwdPc6/nLe0gTIRIxGCJSicypM+giOsO96mkrYZjiugVsAObW6zqD/D6pg6qIjMyKRmRD957ntAhBxnUzOg71LOJ5z7sghEnEU851KzoPxPOejptlcALP+zS3vhVO0ca9720RZkUqrpvACgwPobftQzzHuW5COyHVSQw6VJQSpc4g0dEFmeRO2nTTXKUm1xQawfjOh++//+Hn739eKHCYLzVMJeZfeuCfEtILXstwc02iAfww3zLNeRHXJOpDt9m6O/+7B0bZr4dbpjm3zjWJ+pD/Ld/a4S+GeXM7payH3ZWc5n85qrC5nVyzqI9vrAVWok3gOVf0EuKJDoeZI2bPelxBc/gJ1yQaQLDY3GDwmLczH+RZNy9fu3/ZdLvnvnAHj+d23Dvzm68353hzH27uXD7/sCLOL3pANDmtRZ471p25nXk38LRu8jYvecH55+DgoyPe5ebFJeKD9uYvK/wp6/4pauSXkvDb4R27n9M8L494m8PAE+IWvOl2H10gJ7qPedbn1gqeR1yTaACRyzKevJ0vIG4v3Xme4M8g+HP4eHM+58+Livyc+4JrEvWhLul26N51eGd4c8d9BGTAg5vDO+6gdc7tDs4d0UkZdAeHK+L2sge0rbNijADepfJoDv0ORUM15dQDna36gyrtbhKXPTBzq3Rn89jhmkQDeFJV8s1VM+bsFe9d9sDYqa4+IXMXe5pLa4kdFMVmxSf+3AMnib2p4s75YYYnqB/rHA8EEs8KtebmaZ4XFXp3uAeWd5VV3Lm56T52W4No7Nw5Cl7uBOeCF0Hr8S8XELDBo0pZH+yBsHVWjhIgEmAYtbrBqcOvwZ9uq/WCZw0e0/2VyP38Znpaf+CaRH2o/7PSnVak6nM8rRevre7L4+DFxfExnbE7vAqHbiq5ZlEfb6yVPOee84YhXN1W3mvrEe/4yH0JEnf+6JjWevAb7JSLobkhNdcs6uPPlTR5vONgcP5o3grSb24zCD4FxRcMDjPdLqTzfHlKX/TAKt+rKu7k5YReYVJdtFSt+7r7BxXDzm2tbxxHPbD6Jalc/mgaU93fC0met07zuAfGzsk2uPNzrknUh/rN/HBDqFHtqBcWbWeffPbBZ09+T5cnTz5A5fe06bNykyhnrlLzDdckGoSaKWp1blddZupjDmflt9e8R5BPd/+SZTtgmHoYPCU9sATdDqhnHwbPh+LPh5KfD6W/fTD5OfIweBbzc1LVA5PNO0OuQjyV8imRSNT9k827Qz0LG90UfaL8feZp2BgZEuXgVMkmx2XT45Oy2fEh2fr4umxofBZM0zIwqcA0C6Z1MA2N37Emd5lhGJHOOvM8ddKOQjfN3aqaehZNwCZlrMTtBofrMMzp1AbpEAs8ZVydADAtmpoqXMAi7ThPFVdxqxSxep3ONFf9kHKK1et0uMtPdq9HknE2QqMhhb2ZGWf52dfH6vVI09wprg02LxfkcPw0iGbZ+zIO/akUsfjPBQ7zU8nmdQzt7W/VTWFW0lT1llrW1vz86Z+aQrK56letSJq26qGxR53E71ppalvzc+xRf+cwuNpK6LV1vtLFPDd6h+dDiduWeLa5H+okz9bitp3jSod5tpJiPdQP9U5+Dg4+Gqz+RmFTi2cr0o3d/ExeJavb0Wa1NtEeys+r1cPkWP/go0fgWEDyEDb9yMejsO2/qv3p1vKzrfOy2jwHV5M0pdWxsf7R5NVY/89jybErYJ4Efx6Ojf1ch2ev5Oej5Cq8/+jwanA0uTrWP3aYPFwFmsn+sdGrw6vD0ZF6PFvKTxb9maQzcPVncN3q4aPk74Bn8tHV4aNR8OpYf/247ZH8HDwcQ1kIDjwcRTxXCzyT4M/6PHunvx0bS16t9v+8Ora6Ogo8V5PJJIQxHbVjyeTUPclPaOoodLGD/auHg4fg3cP+w8PDwUFIUzjsH4X92h/uIT00SFNhxALsD9LHcDjYn9uvybNH8rNFPJj5Su/kZ2s8m45bqU5X+EgPzVea5onO8uhBng8lbh9KP/RQxpVu0gmDnUMXrcdfjTaF1eaqj3WNjp8c10knJ3SRSafOORnRTUxKdeOTtEkKpgiYnGBCZlRzWiSTNlaT+ZuttHR8/3/axrJZqKdZ+6qJa6H9H+McXeckYet+47J3FqFQaIz+ygVRZyQy64w4O3/agPxXj5CB5VsWvu4mIvR5bx2/5EE9/p1FWMDWb6wT1dE8Ox1Jr/64JSyFcTHC9sUHQyyc96Y7MAlvYOsrlu9p1Pnz+6Tv/DdZIuw+7uSXVkACNDv6045XOpOB/VdWs1QkknXu5CXJq3/Yq7NEuD5lMUsjneT5zbXxdprg0nc1Pitv7+mVSlHHBJHu/yBkLaXEjCiGTSXco6JbZzBtvjZCvd6hc/Nfvdl1WITG5VKelkXEzVM8jpqyr275/Gx7eXZK90mXjfY/7gpNDqHRz+So3W60ZO12ockktPjBzXa/KRozwhBTpT8yyIZEk9O1Jw2TQ01ho7nqjXXOk38CZeBfy24ZHcLdxWwUaC7uZg8svx1kPR6PxRGNWfyOxeziPvwEpmzlvX+UEenIeu76i49u4LGS8bSWL24GTda2NbL08eYA5aA/ZnGYYn6HULhvEZqWhcaYPSa0xDye62uL5yALroxmaU/7v1dWBikTtxLn+c0ncy0tvKX7aa2C3zkobPXT+dX3jMwDhte7MT9k6DO70AIvWfsB8PV4FncXFz1AX7jF8BSaFr+pGGLoc1gN/7tS7Tlr9BN2O8qTb6s30Zb/lh9MgKdx7TfLmtG0ZgJnmuxr9jWTZ9fj2Vo0eqK7fuPyYqGXsr+7+Xfp+3L8VPGYNZr30kt5OU8FoJHWKwqbehVtdVSM4V1BGdh3aa6etRgyeWIxjyW2HLN4PKbd2LLJ4nB4IKoLQ87im1ykSKXFkKl8bt42w/akjKfL6/W6FC4Xc3B764kMfCZsK34wXLFT4FkzPw0fRUtGRyNTcgMmejHC5tqT34dSUtsSY/40SO51OZOv6wy3hW3zwMLe9srA0vkCw3PgpbKEp8KmD7u0NiKsUIgVxAY4jFDARgwbAlUSMzXhPT2qEGbMsBs+hS39lldRqNaAP185asi8vJcdtwheGFkjhr7c1GIqAnIowpDae7p0vvL23Hy2tG3eO8s9I1BXxvNEoQh7w66wVp9xiWwuvT4j1mq9Ln5GryXCeq8NeUyvtxGnJ3oX+BPMYUKrP8lM2Wx6LwExILKFvXpbiUNr5afUcxuDUtSQgpbvJid1+atup2cj9LNXB56eLzzdO1sZeIqeKZzzp/kmT75Y77K5bHwXIePbCIU+fEqItS6twgU7Ci+MKl4w2sB46sqEtYqw1mZTADdigyAg0MHRevShEn/W6G9/XWyEZyW7gja0fzs+Pq7K8VRNqHI8F1bMK8BzuyZP4tRlI8B1wNMFvgxDiHpPIBUVU1rtKcEnNsSKTEYLRpvNtqHVe7VhBR/i1ubVgj8JL0SuwuYq4VkjP9WGWP24rUA0Frumd6ArUufidmpqwgBd6t+ZGN1e2qZ5Lj3dM98atzatyxYWg7tOCVkYXAM8T20ZRTh8CnkITOBFa9Pz+V5XxpVRiMMZ1HF5xWGFDXKT8Cq8fIW21J+1+9tXi6h3sYDU2wIv+f1CiwkJPT/Tr/rtFpB7FiO8Z0J2i2XLJFw2CR3Io/6vGF00JJodz32JKkfq3DywtzRgBomwZ2ZMfzWU8UQOIiBuTzMywps5yZycuvQK4pQ4zXjDGf2JHiWcNqMnNrRaWzjDBzPEc+YUivZEj/pf8DJUEzeWnwjOZ8BtP7u2uxwTLmaXo9Go0OM/yGaRy66Xs/uWtbXo8oHD7we75xpGUPRD7MOPcyDN/WWDsqAZlAtVVAI9rkjKxhUxQUAPKoa+1IUO+C4x2vKh2S4Ceig6HOE9hYKAPjkMgwsyi11iBaobprtkAr3f8LhC++NLkDnXUeEykjxr9qwwZt+3+59BHq4JjQ67Qxi9FtrX/NdCT3TRgyoaQeVuiarOW1TV5BDohJ/67qqHoF+2NfDBBnj2GWCecu0XZj0xxHPZnrU8290F5YOk3zJI3EXguc/wtMNGGD0QZoeqp4P838zVvPmp8s48+XyiMT3UwD9YDN98mwWefofRsm/y73uEIP0OQBzE7CD9YsKtXaFn2b9r3PUsRoXZLZiF/nj7yqrypXmlHOZ/zc1X7ubPwqY2GuGpi/R9vSjcNXocax6hHSI2CkmJehrH7hpy6m7MYTfFQPctrjlgirovqjWlVupkExPTsxOT0+sTQ9NDE3pdrgXN8AwTeZoZ1NfYalTN1WtoXgYY+t5ycwAxXZssZcsLwgMYhq7veufZOjyRss8VhULryh+GxbDVi+vp/oZ59vX9UKEaoo618rXca5P9u9vWTVrkSXhB7Wm9epB+oARGvK5TkA5eLf9E7DoFMej1nohrfbyhuM1B+qzqCnUJjNctLE7V5nniEmeQtLOdnIhPwqAJT4iMC6zaMGgDPRKERG1/NrG6LJV+WZOm/9mdnVmPJ2gEPj/jBe1+4gIdC1rJxdeKFS5QTSB/tGLtqa2mP5uIWwT5v2zdytLyZWv3hK7JU+wFUuDTsIvmqQ8jngS4l+YJkxZCX2Oe2ti6SSkMzmxFf8Rg66sWV+Nrx61Lm4EQtXnDiKcro7eBnNVnvISWyJxoT8KnNv3NqXU5mshPBurH1Vxqyrb88OPaPBVhGxQQdTCkEITYRRDokABJD1M1eLHVlgtN5WcOXzsqXBr9U+ur4/XGTwW/QhAoqu5W53mHU7oMP0TLWNq/a8fd2ju8rnmn/+a9WitOTI1RZ1v+T8bxumZ1GCbzq532tTb9Q4fzderqkH+FVIPxP8bb9X8rrUIMQyBdci+lhV9h5jdcE5W75GeO6I8HRv+79p3zMa1//JFeC+Xk8Yn+o8daKHr9Y1QYc5np8anox8Zq5v7mRCtNc7bzqRgGg1wORWKQ0CV3yBRJzlw0RUQ6eWM1mb/ZUtNa+pFaBIu3q+mb4urpmgblumhdyda3Szl7xPy0iM0bvslHaq4cdBDMfe1YunugZIPFW66VwzDC4v3eDPBVG+x8VQXYvK8dOpew5VnKHYFOTGXrzs1S+p6s7HxXBTp63lsZpCIueU6Lptjqhuh7XG5w1OMaRCqWvmmInXO3b4FBxNZzEuUStcog5+rOdnI2e8C23j9B0hxmDU1Vbyns2nod3b//c1MINFc91TX3w9DMYJ3DzMetJFhbr6PTYILOAQt0zXWRZTzBBdVaW+OoHs+uuc61hCeGk5SvSmt9ZdSq1RDgeHWeH3fNdcslPPG4JpCgkE8LfkVJhmsw+pjZ4gkBlnM8VtiZCZC5fawsJLBuzE8sRM1geFrgo3A8QIKLMBz3kQEfeNBHUTjme4/0gc8T9H6AEuDMDgkuptLAk973+QIlDm8tbjuVn3EcPIP5XoR8cSoUClAzITKdIlO+uO8FFQhhcTIFvsRFVCgNJg0VCpApQYIKBHwa8gWJYWBMUek4mSjGeYv52SGeCRwLxRN4YoZKz+AJ4Bkg0/hMKJTyJWYw4DXjiyN/CmZS4OhQIJAiBVTKRyU0JIpbDD4g0KTQr1Pk2Y1xK0gjR2iAZ6AKz3iohCcVRw4m05pQiiR9qTxPDPHEynh2YT+E+RKkj4rjL2bwOA4ZB9FJpkN4nIzTPPEEGaLjFrKRCuCaUJrCE74Eys0QngKepAYPBdLlPLtzXPGF4gHogzCM0oRwQUBD4emQBoITTAIKI0MBDfS3gXQax0MhisTTKR9Tk0oH0LhCpgJQFyOLY0yX6oTcWJHTAtCDQhTndQEML3iIwgo6AaNrF2sKqimI7szPKiBL2kyFqFpVqwH7uGvun1CNZ8FlZS7CBFhFlXo8uzI/C63LjfS+W8ngOfVXVQ6X8eyeedkNHY8Qwmm9l84d5jKRVnVMrZAvTVehfFhR7BVkYMnf69L8JMkATgXAWSSj90iKAsmHspQEjQ/9Lj1kUFQaB4UH0i+eRi8YqEN4AwYbjAyQ5Ty705+JdCgBwwmMnGkqhacELwLpBAXDJy5I+URUXEPBHgZjygswpXxpMgUSEdeQoTgMuwENmYKhVlMa7d2anwkcj8/4kH4VzMRh3E9heGiGYniCTArMgGQCxSDQgCkBrguRSMMHQuQMcAbdDz9PKF02X+nO/hYRYnhiMwkKeM4Az/cYnvEZKs8TS4MJpxKBkA9Jfpon+B/DUmmKIst4dmd+Fv1ZwjMAXAIvijzTFEnHLSijAECgCdE80xDDoTSjJYpx25V6SAAaLoShnigEaegL0Qeg40gQ8TjodhKjQNPhaVC9IQEVB/2nwTUgdkHsIRkYJwWBeKB8XOnO/MTy2i+/OsCog/zwktd3WNGYe0uA5QeeG+NKF85XOoBuzc+28+zKdZP7zbOT6/FdFLeJ95pCqLnqmq7ph5xOpW4iopRCiUxIoeiUzgkw5c15k06JTCqRStdYTeZvtuQR7p9zyg44e84py8+35/A5bZINFnly9pxTlp9vz2F+ymVs+pO75yg+lPyUs3lCD3fPrX0w+flQ+tsHk58qNvOTOz2knmXxvrtcPod4RCpRGuQGpYQuSqVcDoeMKW+WI7OkYL5zTUmb7z3YFM8N0fSITCUbmVZtQBkZUalGNlSMKW9WIfN0wXz3mu08+fb/AZJr89OcuOZgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(filename='data/index.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add required packages. \n",
    "import numpy as np # Fast array-processing package\n",
    "import pandas as pd #  To deal with structured data e.g (tabular, multidimensional, potentially heterogeneous) \n",
    "import matplotlib.pyplot as plt# For visualization \n",
    "import seaborn as sb # For statistical data visualization\n",
    "from scipy.io import loadmat # library used for scientific computing and technical computing\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we intialize the centroids for the clusters. write a line of code to **randomly initialize** the centroid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_centroids(X, k):\n",
    "    m, n = X.shape\n",
    "    centroids = np.zeros((k, n))\n",
    "    idx = np.random.randint(0, m, k)\n",
    "    \n",
    "    for i in range(k):\n",
    "        centroids[i,:] = X[idx[i],:]\n",
    "    \n",
    "    return centroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the function to find the closest centroids. The input to this function is the dataset X and the centroids. Write a line of  code in order to **compute the distance** between the data point and each centroid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_closest_centroids(X, centroids):\n",
    "    m = X.shape[0]\n",
    "    k = centroids.shape[0]\n",
    "    idx = np.zeros(m)\n",
    "    \n",
    "    for i in range(m):\n",
    "        min_dist = 1000000\n",
    "        for j in range(k):\n",
    "            dist = np.sum((X[i,:] - centroids[j,:]) ** 2)\n",
    "            #dist = np.sqrt(np.sum((X[i,:]-centroids[j,:])**2))\n",
    "            if dist < min_dist:\n",
    "                min_dist = dist\n",
    "                idx[i] = j\n",
    "    \n",
    "    return idx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this function we compute new centroids for from the clusters. Write a line of code to **compute the centroid**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_centroids(X, idx, k):\n",
    "    m, n = X.shape\n",
    "    centroids = np.zeros((k, n))\n",
    "# compute the new centroids   \n",
    "    for i in range(k):\n",
    "        indices = np.where(idx == i)\n",
    "        centroids[i,:] = np.average(X[indices],axis=0)\n",
    "    \n",
    "    return centroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the function where we run kmeans and optimize the clusters. Call the function we defined above in proper sequence following the flow chart above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_k_means(X, initial_centroids, max_iters):\n",
    "    m, n = X.shape\n",
    "    k = initial_centroids.shape[0]\n",
    "    idx = np.zeros(m)\n",
    "    # 1. Call function to initial_centroids\n",
    "    centroids = init_centroids(X,k)\n",
    "    for i in range(max_iters):\n",
    "        # 2. Form intial cluster by finding the distance from centroids\n",
    "        idx = find_closest_centroids(X,centroids)\n",
    "        # 3. Update the centroids\n",
    "        centroids = compute_centroids(X,idx,k)\n",
    "    \n",
    "    return idx, centroids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "128 128\n"
     ]
    }
   ],
   "source": [
    "#from IPython.display import Image\n",
    "from PIL import Image\n",
    "\n",
    "#Image(filename='data/bird_small.png')\n",
    "#w,h=Image.size\n",
    "im = Image.open('data/bird_small.png')\n",
    "width, height = im.size\n",
    "print(width,height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "image_data = loadmat('data/bird_small.mat')\n",
    "image_data\n",
    "A = image_data['A']\n",
    "A.shape\n",
    "# normalize value ranges\n",
    "A = A / 255.\n",
    "# reshape the array\n",
    "X = np.reshape(A, (A.shape[0] * A.shape[1], A.shape[2]))\n",
    "X.shape\n",
    "# randomly initialize the centroids using the function defined above for 16 clusters\n",
    "initial_centroids = init_centroids(X,k=16)\n",
    "# run the algorithm for 10 iterations using the function above\n",
    "# Bonus run for 5, 10 , 15, 20 iterations\n",
    "idx, centroids = run_k_means(X,initial_centroids,max_iters=10)\n",
    "# get the closest centroids one last time using the function to find the closest centroid\n",
    "idx = find_closest_centroids(X,centroids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "49152"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD8CAYAAAB+fLH0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXuQFdd9579nNKBBgHkEGA1vAgKErAGrNAgIcrkkO5FlESmK1yVv5JBYKdVWORvH66pESv6wtyqpiivZON4tr7OUbEdZeS17Zb2QFdmyZK9MMYJBChos0ExmJAYYhmEUIWxegmHO/nH61/3r06dPn37cO3dG51NFcadv3+6entvnfH+/83sIKSU8Ho+HaBrvC/B4PI2FHxQ8Hk8MPyh4PJ4YflDweDwx/KDg8Xhi+EHB4/HE8IOCx+OJUbNBQQhxmxCiRwjRJ4R4oFbn8Xg81SJqEbwkhLgCQC+AjwE4BqALwKellAcrP5nH46mU5hoddyOAPinlmwAghHgUwJ0AjIPC3FnT5OLWWTW6lImHfO88AEBcOS31vUbi3KXLsZ+bps4IX7dIdb0XRPJ3IaZNzX9O/R7x+6LfN/neeeO91Dl/ETh35gwA4OIlNVmevTQWvn/58uXEZy6PjSW2NSqXRy+/LaWcn7VfrQaFRQCOsp+PAbiJ7yCEuB/A/QCwaMFM/PDrn6nRpUw8LvQdAAC0rLo+9b1GonvodOznliVbwterR18HAPQ2X5f6+fbl+a1Y/R7x+6Lftwt9B4z3Uqf78Bhe7dwNADgypB72vYNnw/ff/dUvE585e+ZsYluj8u9vvzPgsl+tBoVMpJQ7AOwAgPbVV0/IBAzbw1sG2/H4ey4DBD2w7W1mJaY/0Gn76fCHvGVJcD1Hd8f+b1myxToYhNdweCwxMHQfzpiBg+O2O12t299KXUMwoAWDw95B9aNpQJis1GpQGASwhP28ONg2qah6MLCRVyHoDzv/OevBTxsouodOhyqAFAAQDRD0Hh8cuGqwnjMYBPKqhnDwaL4uvKYygzWd/wLp3H3qv9kzPxAODLNnfiDxuYmkGLKo1epDF4BrhBArhBBTAdwD4Okancvj8VRITZSClHJUCPHHAH4E4AoA35JSvp7xsQmNPpPbZilXG9e0Pz+PPqO7Yvuc6T2+jVRAN3u/ZQlS4aaE07UFM//81rUAgJHhN5w+Z6Nl1fW5lRZd78ZFuwAo3wIphLlz5oX7TUazomY+BSnlswCerdXxPR5PbRg3R+NkQ5/JbWrAdeaiz6cdh9v64wmpgSrhCsFVNehO1bx+itb1t4ev24P57EKH+mzfuXl459TbAIBVV0XLn286HXli4cOcPR5PjEmtFGq1ZGjD5VxZKsG0zm6DryaMt2pIo6f/Uvh6zcopsW30cxqkEPIrBuXG6j58HVY7XOPwa8/G1AIAbLx5KwDglmO78GK/8iUsbXsPQBTLMNmY1INCuEbO1ryLBMpw+Jo6yVO+PFcFJqclbcuKO5gIuA4GOjQY7B3dho3NOxPv60ui9Pcv8vehwWH4NWVGrFsMHBlSZsP0BUEc3tCecP/pM6bHPj+Rlyi9+eDxeGJMOqVgio6r6rgEzdokSVtWXZ+IwLPNTtw84M6wNDPhQt+BhFnQy4J1wmM5mg4Hj6n/1y2uZj99f2L6gg7jfrpC4KaFDfrcrIHHAQcTJFQMPOw6w6QkZWA2I9TyZMsSda6zJ5swezAZyERMVLXglYLH44kx6ZSCDa4ghs6reb5tWq9xX9vMn4jnN8Tpp6kBnSg0N9pmsoX15ccLR3cDBf0KfObXZ/d1i5Pb+M/0WX2fPNDs3nVc3Y8TA92Jfa5etiF83bFwLPY5/bWO7lsAivt9uG+BnI709zwyNIaNi5QvgSdOTXQm/KCgP3A20yGPaaE7E/M6rboPj4XeciBfVJ7zOWqw0pD1sLsMBmlmA6AGghMD+xPbae3/zPzN2ScAcHrZ3QACU6IA3Tmdz63rbw9NC2JpW1O4AkFRjhTLMJHx5oPH44kx4ZWCbZTXZ3sTZEZwuEmRPP71obOKZuq0uH59Pd14HcwUcVEIjRKTQGrg7Mkup/0jUyFSCVwdXLuQ7o3b2r+uEPaObgu3kdORmxG9jrkXRJrDUWdpm/q99u6b+AqB8ErB4/HEmPBKQcc087cvVzN/9+ExloWX/Gya05HTfXgsnN0pO3D16OtOxUQIU3RfPTD5A8o6Dqcv6LCqBZNCICfimeDnjoWROiBfAfGJRf9mPO6u3cnl26L30jVHglRDmBdxFHiua/JFNU66QcFGfPUh+X7WikQaeQYEQH15de95T/8lYGXxCLyy2FYkXNEdjNypyE0FPggAyYGA09R6o3H71sAaoMHBFrvgmradhe5obFmyBddeVOfqOxdfQZnIKdXefPB4PDEmnVJwneX5frrJYXM+crVhc2BmmRS61LVJ31o6FMvEGxAmhQAok0FfarSZCpzIbDArBVIQXDHsHd0GAJgF5XAsmmNhgqsEbkZcOKpMp1tWqvde7L8SQLxk20SLbPRKwePxxJh0SqEIJnWhq4AR9loPlDIpBp6bYCuBbkIPlGpvm1VYLfD8haJORQrQoeU3HZMzEVB+hGvXU2V/N4VA/HDwGgDANoNDmEOK4cO/cyN2vlq9ouJLk7pPgXPDZvU3fbH/FQDep+DxeCYRXimkYOtDoCsD01JW9+GxMP+hHfFgpwtHd2PgorLDyd7l5dD14qiutRNsCuDgsUgp5FUMJoVA6gCIFAKF+FLIr/IjuK80mOCz/7Yb7PeB3h9bpO67admyKGkqgf5m9B24ZaUqwPLmkcpOXXf8oOAINxVcYuX5PqbBYdlU5aCi/gLtbbOQ5iLtHjqdGBhM5oQpmakKsswHggYDikNwdSpmQQ97lnmQNmj09F9ycjba/rbcfOADhL6/Kc2cCrBMFIejNx88Hk8MrxRgX1rUZ4L25U2FuxmFxzA4DtXPuxP7VUlZFXFw9KbENlOkYlUKgcjrQNSXK117R9j+nrx+o83hSOnVf3RsFx766cRsiuaVgsfjieGVAuxLizYVkTcn39TlyZXxzIgkX8L0BWMxByOglh31wKQq1EEWvZ1PAgBWb74rdZ8wRLqCLlOAXSHo/MGnt+Kxff8CYOItTxYeFIQQSwD8M4BWABLADinl14QQcwF8D8ByAIcBfEpKear8pVaDTfqnrSK4oD/oPJrRdFy9UWtWuzYX+ApDFQ5HcjCGhU+OR++R2XD1smQuQxYuDzTtY2L15rvCz/L9diJ+PF7xOSp4o6iiHR2QrPrM+eSNMwEAD/00GhQmgtOxjPkwCuCLUsp1ADYB+JwQYh2ABwC8IKW8BsALwc8ej2eCUFgpSCmHAAwFr38lhDgEYBGAOwF8JNjtYQA/A/Dnpa6yQlydg5T/QGnX9HPaDKNHIcYzHeOVg1tWXR/WdYwUhb31Wt4KzK7oDU14iTEdk1MRiMyFcNY+Hs3ephndRNr7Jwb2x85l25e2hzUdg+1btyQrN89vXRv+LfM6jLMKrxD3dChH8WP7fhVua2SFQFTiaBRCLAfwIQB7ALQGAwYAnIAyL0yfuV8IsU8Ise+d0+equAyPx1MBpR2NQogZAH4A4E+llL8UQoTvSSmlEEKaPiel3AFgBwC0r77auE8joNdYGLHtnIKtfZ2pOauuCrivoGr0gCSuEsiXYCq0SqzefFdi5ua+AptCyPIb6Mcyfc70PnliTMFMY8P7wtd5FAJXB66l2ujv/ckb1RW92D8Pbx5RLWkbWTGUGhSEEFOgBoTvSCmpaN6wEKJNSjkkhGgDcDLrOJfGWhLpyi4p0EPnV+cuiFIWHqeQF+6MpN+WQplblmwJB4gqnYVlMA0GJOVNZoG+jTsEqyB8yG/Ynohd6O18svC5bMV1TIOBCZsDm5KljgztCrtUN7LDsbD5IJQk+CaAQ1LKv2dvPQ1ge/B6O4Cnil+ex+OpN2WUwm8A+AyAA0IImlL+AsDfAPi+EOI+AAMAPlXk4C7KoQqVkKY2aBtdR1YUY9kSaiYzguPiaCyS+5CnczJ3+L306JfD1x++58ux/fiMnZXEVBT9uHw5sug5TcV1TA4xF/Vwoe9AlPMSfGeeBPCba9sAAE/s6zN+rhEos/qwC4BIefvWosf1eDzjy4SJaMzqz8D3sykIk/1oUyXh6yD4Zeh8tExJPoK8hVs5ejs4E1wBmFq42cijBIi+c9MAILEUyNHVgYlaqYSXnvh2+JpHT+Y9n+k7RWzYtApA3KnoshRJitJ05Kws00ZhYlylx+OpGxNGKZgwqYIsP4Prqkba/qbPm/wJphoIeclahbCpjCIKAVDLkFdb3rd5+F1napesx6wekXp+heu5KR+CL01ySCEQroFKev6Mafn5hs1bwmKy02ckMygbZSViQg8KVWAzI9K2EWUSnMqybnFyMCizhGlr7Fr1QGB64PWGMvR4HBkasyZc8fPTOagSdFrPCHovbWBwZf/LfRgZTh989fgUNXCoYi+/vvTXAQAHDh4Ilyfpf2K8BglvPng8nhgNoRRGRy+Eceh6NlsWtpmcVACfyU2yLo/z0bkcm6PpUKZKs83RSE4tWyk1PgtzigYBmcwCUgXNTAm4zH/US+LaBYBe53HWwOOhWqDZvqn1xlA1tK9QJkD3W+nHT1MJ+1/uC49n249I6y3RfXgM7YbvGkH9MN7U1AGHK4d6qgavFDweT4yGUArTpkbhw2mKYWT4DauKMGUvUp4CzfkmlQAkA5X4NiJ6L3mellXXj4tfIQ8mx2PfuWmhUzGvOjDZ8gSpg0Ov7QlLMHClQgqH97+wFUzVUbN5tt+gfcW7AIDut2an7pMFVwz0HaNOVAAwC/n+7vR7vhrcj9mDUScpKsYy3g7HhhgUKPeBd4I2PuSWwhg8IpC+bDSItExzE0R8cNBTpSk2oa1gXcYyZEUq8vdcVh14HEKZwUCnuetBAJF5sLStKRwAyJxKG5jT4LURCb5CsP/lSN7TA2wbBCJnZNJEyHJM7g0a1BC2FZILR3fjwmj8d77QdyAR0zJ3zrzQlHhin0ouNpkNPFei1nkT3nzweDwxGkIp5EVv4w4Aa1ZGUlQ3M7gjMW9ORZFU6bLUqqAKYYtUzMtLT3w75kQEkFAHHFcziysKPb+AKweTaghbzhlUTbQtUgW0PwbzOX1PL7s7VvLNBYpp6QkcqauG9kQ1MA1OR32baZ+q8UrB4/HEaFilQLM92fQ9/Zdw6LU9sX2WtjWFueo2TGqAjmtqHGu6DqT2b6qOPDP/kaExa4EUE7T8WKbGAdnkVDNhXXNXqjIostxKxzApClIPprbwtms1KYa0qEpbwRgTre33Aoj8EtzvFZboY7+Lfk+434UyKH/8xhDGk4YbFPggAAA7O7vD91Zdpf6nh8FlQAD0UuzpFZRM8Qqu6dl6jcZ6kCeU+cz8zaULnvR2PokZI50AgHXB32Dd4sixS3eKnL7ULzMNUyu3bi2FnA80pngTPkC0TaNXcYcg4BZazQu1RFGR6Y5VfkxyUrYG/w8b9m9ZdT0wtCu2jZuIS7UVifEqDe/NB4/HE6MhlAJFND7y88uJEmDUxXf6go6wKStf3zZhW7o0lVIjs2HofPa1ppVisykEkox81nOR1q4ORxtVmAwkqWeMdLLGMEoFtCxJzvb099maMzoVQJiiTvCGIfR3XT36esK8sC11ZqmEWBXoYJlx14B6b83K9EjbTywCyGGpmyqt7fdiuPsRdXxahmTfHcr1OIhI8U6/qNTx3DlKHXul4PF4GoKGUAon/v0s/u6RVwAA2za3A4jszbjD0c2HoKM7LdOoorybiwLI2qeKys22rMcsuDIAgBnB9qVtTaFC0P8+aegBQRRlmAcKRiJ7/dTwPuPfKix606lm3Cx1pDsVu443YV2zcmbT7E2/nymwaWx4X+hg3HZDUjGQE5IUgwm6n0WwBTSVyZvwSsHj8cRoCKXwgSslbln5Hm7YvCW159+Fo7sT4csmbGqALzly30De1vJ6QQ3XFYeiCsF1qbLv3LREYJKtB0NaTwVSCMS166MW9KZuS3phkjj5lYFOQl2sWIXut5JZjCPN6u/RsVApBdvdNvWcWNe8J/wbmBSC+Tpma9ehVj52vno69C9wpar/LU0rMGEGpeX6sygTAt0Qg8JVM2Yklhf1pcllU5PvuaZZhwNF69pQdpo6TZsGB30bf69Wy482B2NaYlMaWY1T9P1MzkSbqWAfEGpHmOwUOPq4RKeBi5YTsxyNVLzlLOtQEn5n2CTzfDditLbfm7iOWbuj+pE7odK7ybQwdb9uX94Umj0XpqrrPGS9WvMDX2UehDcfPB5PjIZQChxd/tuWIU3p1GkmQnjsxP69if157Udbe/r08i71QVcIeXMaejufDJeASbJeu/6mhlQGadBMvX94LY7u+xYAoClw8NF7Y8N9UX6DAVoe5MqMFGrXcfX3PzGwP7xHt3WobQNHd2NWsMzc/ruqh/JL7LiURUmKYdZAlLPDz6U3GeZBTOOxLOmVgsfjiVFFg9krAOwDMCilvEMIsQLAowB+DcArAD4jpbxoO8aF92Qi85FmK1ugkmkm41mQpAK4YrD5IyJVYF+arLUvQYf7EXiZtXdOvQ1A5eRnkeZw5AoBUPddvzdVqwNbT0Yd12rKsc/88v8Er/J/lgiDi/arNvLv/uqXzPGnchRu62gKnce9Tz0T+/yalVPCYizR/W5Cx8Jk0Ffo34I5iInOr2Oqq1BFrYUqzIfPQ/lGqITMVwB8VUr5qBDiHwHcB+AbtgO8e+ZcmONAErjr+GUAwL03X5H7gkzJTmn7jEfRFBOmJi98m15zce/g2XAwsJkNttWHGSOdsZUFouggkOdhr+qYLoNGU+uN2BYU8OFORzKd1q1MfsYlr4TXyaR8D3JW9qADWBbff/Xmu7AmSLWOJhZWtUsrwKKieWcCAB4LFjcGh4as6dNVOBxLPRFCiMUAPgHgoeBnAeAWAI8FuzwMwB5B4vF4GoqySuEfAPwZaDhTJsO7UsrR4OdjABblOSCN3i5OM+6UdIlatDkhqyJvfgMnrfnL0rammEIgXBUCEFcJFIdgcirmVQlF1EGeepZZ5dvo/G3TgKMp77XH1ISKK/hvnfz+xVPyAV5TUs3KeweT761bnGzIQ3+naxfwcm3R946iHHuDJdR2y+/Gox3JlBgcGqp5DccyrejvAHBSSvlKwc/fL4TYJ4TYd+lSspKSx+MZH8q2ov9tIcTtAFqgfApfAzBbCNEcqIXFAJL9sQBIKXcA2AEAM2fOlLQ9yn3I50vIymsA4j0bXKMXaw1XBXmKrKQ5F10VAhB3KtZrqTFv1eu0/ct251q9+a7w3tDsvm5xsuI0zZtL22Ym3jOVmzP5aEj9fvFz243XEqqhQL2Seuvpv8Q6Z10JQDkSTcVcq6TwkyGlfFBKuVhKuRzAPQBelFL+HoCfAvhksNt2AE+VvkqPx1M3ahG89OcAHhVC/BWAfwXwzawPXL48indOvY11Gz6KrVviHlibAujpv2SMHa8lppmJ+w7yNpWlWee5Lnt5NfIl2JYfXVcauB+hrEIwNfklisziJh+C6ThFFUKUvzAbLw2o+9sFlVH6XFdnGJhEmELNTX4jUnnTF6j/Ty+7Gy89+mUAZt+PrcyfqTAxKeidUP0ngdqVeK9kUJBS/gzAz4LXbwLYmOfz06c0YeOi6bgr5/JjkQFBdzCmmRF6c9CswSDPe4DZmahD2x7b9yvrsbIGAyDuVKTfuXW9fUAgR13Wsl9a6z5KUEpcryXOo4qmOnS9dP2m3hHtK97FcHA/aPn7zPzNeK5L3S/uTNTRBwIgcgpS3sVXH0vmnPDl0E84ut8PjqqBvCNonXf1sg1484iKmOBmhM998Hg8NaMhch9arpqBa9ffhJ7+S+jpLzdTVGVO6DLWtTBr3mXIrAawgIpmo9blNkxORZLD7W29aFmufqe8UYI043JFQGaHyfzgRVH0Bq0jw28kuiSZoPtsup9py77h8iClM7NCvKbeER+7847gJxWN+MjPLycchc91JZcrCV50hvjqY+q6TwzsN5oNlE49xiq7knqNupxFxwzT2BfeFP+5hnil4PF4YjSEUjh7McpGKwrlxOfF1N/P5GfgCqFo+3hX9ECl2TM/kMhz4LOQXSFE/QxdFQK1Y7flgdI+JqXAi5F0I17GrJUVLTG1eSenZffheOYgp3vodPh78eAhUnenguPSte1/OVlqj/sZSDH09D8Rfg/p+3TX3VvDz5gcgATPpgSAd069jauXpe4e/S6WILqzWuctQPmG9g6+GN+vYodjQwwKjQaPZxgPqBnI7JlRR2KXpKcZI53GwQCwdW9W8ApDLp2cXbHuv8Lk6FTb5p8Prq11baz6FmEamMPfMTAVqDpTU+uNGArkuq0O59Yt1wO7lflKDzk5IU2TTtfxptggwFm34aOl+2wAZrOSvguDQ+kNacvgzQePxxPjfa8U1IyjZiBTmnYo72rcAYrHJJBCoHRZrhhMzituMrjESZBCsCkCTpEKzGWJmyX0Wsn85596JqYa0jCZMd3dUe8IvYvThk23R/fkCVVWLVIM0fyp9yYBkmadSSXw9nWm+Bv6/q2hDSu3JEyWrVuuDzOKa9Vs1isFj8cTY8IrBZOtp88iNHtmLYXR51qWbElUaqbPVuVk1PP1eSakC7w8GNmdXCWYfAm0VJelEMZDGeThY3fegeeD4Hn+N7NlbEYFVqMsxdVagNrwa8+GGZVjVLU68DHs7Oy2FrWxKQSCV3guA52Lgpg4unoo4mNoiEHh0sVzuVKmTfD4hIQZEDzYq0dfd1ojB6IvW1T7Jluu5sHkQKLVBr3KDv8S8pqKtqg7ggYCDvf60wDR6AOBDoUJHw3+Tupn5US0rbLwmo69gf3A05dpYNmwKX4MZT7EzQbTd5W3oDu97O7Ye9tumGVcccnDyPAb4UR4cGZtGtF688Hj8cRoCKXAyaMYOhaO5Y5gdImU09uhV40pRfrI0FhixOcORoKbDCaFcGrR76oX1Cz3fPrMpJx51SsELuOL1FisBxs2rcLzTylnn60ZEDk8vwDgh4PKNDA1ktHpOt4EHGeNawG8NGCPuCWnom2f+a1rMZ9Ky3XWphGtVwoejydGQygFSp022c4m9Ca0aZRZPqyiDbwLpjJrJkgh2Ohtvi7sQkSzDm/zVotCKvtf7rMGBLmUaxs6v9rp2mgptW1aL9qClhdUgq1Ic2C9hF/34bHwO2Ncrnwr6DhlKTvqoiJ09GVHU+Qkfdd5rxObw5FDzsd/f/sdp+vxSsHj8cRoCKVA6GpBR58tbbHoQNSbzxTQkxXkk1ZElb9XJTa78J1TbwNXJQNV9GujfgFANLPQDNqotj1AWYzZM30ba4hFdRfo73ih70BsaVHH9Pvz3AgA2LX7QPSdQXy5EogXaNGhWgl8STIq3BqR9Z1Ng39u76jqdjVj5CF1NTOr7STVUIMCkIwhJ+bOmRe2STvjmDx1NvjtBi5GiT3Uho5DS5imKDl9ADh4LPkw1trE2LgoGhD4UqYp7fo/3Hpl7LO1HgyKyHa9kAp/8IoWWXGt+kzw+8JNl13BV+BgV9BKbmgXAFU1iD4Txjyw5V49/iCrqW0ZaLDpna8qRmHwJ6n7Fol69OaDx+OJ0XBKwQWTE9K0hEmlrHA82taFeBGNjoVjQCjNVHGLsye76uZgrOLzdK0bb45UQiObC/qsXkQd0JLyxpuj1GabAtExlWgzKQagKzw2OR+5YiC1oCsDbjrwFGi9KMualVMKmxTE3DnzKjUfvFLweDwxGlYpuNQP4BQNkzYXd7kJBw/Ht6xrTi/LlRdXlUC+BJ4Xwf0IkUKIZst6KQTbUiOfsU2ztV5YtQqfQuv62zGgfTbrOkyFafWl0VAxDHWFzkeuGKKwZeX8S1MIadsuTAWWTVWvBy52JPYfDxp2UEhzOHJMA0dabbyyhKYIolRlU/5CXsh5mhfesoyop8lA5+KFWoi24Lm60HfA+JDrMQAc1yYv5Bzmv7PtszbTwjY4UAzDq50AEF+ZUL+HcmKbVhpc4E1qyQnuOjhQpOShkfN4p8I8CG8+eDyeGA2lFEyjnCn+n0hXE+nRkK7Y1MaZYCnoYNBG9+DhpHlhclTypUybQqDfmbef11VJe9ssa6m1erFh06qEWqCqz22r3Bu56PI+SzGY8hToPuhmhO08HJPzkaIHpy+4hIPH1Eze3ha9n9as6NBrtqa1cZKxMOmOSU4U0Zus21gGrxQ8Hk+MUkpBCDEbwEMAPghAAvgsgB4A3wOwHMBhAJ+SUp6yHefipUthEUoXbOoBcPNHZJOtNria4D4HAMCxaKYwqQbbNfJgJSA+w9zTES/IWgZTrYUidRXI/g77LQTOt+7DY85l7PI6GG3qqEzzWd2/wH0Lr55U+9ByaDsOhPWue4JX5LieYTh2moPZ5pvSHZNcOfClTPrO7A3aOZfxLZQ1H74G4Dkp5SeFEFMBXAXgLwC8IKX8GyHEAwAegOovWQheOYais0y1C6vGVmUneqDT4yUOjt4UmhSu3aR1s4FTNm7CNAAA0QPsWq/RRjighFWaI9OiZZr6napoC5c1IIZO0MCMmTP4g/A9vR1gGvrgoMyV9MhXnb5z08KwfDIV05LabA2BdEwxD6ZVjjIUNh+EELMAfBhBA1kp5UUp5bsA7gTwcLDbw4Alpczj8TQcZZTCCgAjAL4thFgP4BUAnwfQKqUkW+AEgNZylxjBG2oCSjHUUi2Y4LKfvyZFEYu2XKZMCt0JubStCT9+Iy7vzp45i99c2wYTfPmRz3AuDkbTkiFHVwhZqdB5zg3wdf/g//W3h7NwmZiEPPQ2X5cwX3gClQ1TPAYtGbazjtE0W88YIRMhciS7pL0XoWqFQJRxNDYDuAHAN6SUHwJwFspUCJFSSihfQwIhxP1CiH1CiH1yzLiLx+MZB8oohWMAjkkpaRp8DGpQGBZCtEkph4QQbQBOmj4spdwBYAcANE9pzjUqcMVg8y+4lDfLQvctzJ0zz+gktPkgbDYj/S6L2tpSbcqBix2xWQmIz5akBvjyGM2MbZaIQoXmWFyxCjSr67Mkn9kHDNGCrjO3qTa+AAAefklEQVR4uB/b30U95HWu8tRoKtLKFYOrf4GIlkHVUmD30ah9Hfl8XuyP/Ae6LyGPnwFwj3ytIoiOU/hoUsoTAI4KIah3xa0ADgJ4GsD2YNt2AE+VukKPx1NXyq4+/GcA3wlWHt4E8IdQA833hRD3ARgA8KmS50jl7JmziXzxvEsxrn4J12VOm2KwoS9DApF3ec3KKYk28tz2p+Ijbcub2Gx9B8qSmPmZP4DgocwDhpm3qHoY+MFXwmOVXbFoar0R8+mHweTSqE0x8PfoNS1Ddif2Zt8TVhSHVME7p95GH9T3wqQeimbO8s9VsTRZalCQUu4HYFrLurXMcctgq0KTtb1qpyU5HW9ZGRdkR4bGwnOZzmmKYqMowe6nngEQyOBxiGjUzzWMpOQ3mRncxDDlGugs+121il1FZej2Fe9if2A+9FriJmyDT9p7eZoDxUzPYNAwDRTjjY9o9Hg8MRoq96Ge1Hopk9ebfLH/ysT7c+fEcx/6zgF9/er1qiHlu6UcC+ASevrVTEWZdC0sXXo8aV1/eyj5bc5CbmLYaimmfS6NovUYe4erbRZMAWqR2Xje6FTcG8T2RrkvZ40mx3jilYLH44kxIZVCVjFKWyv3RoZCpM9o23khGCon99jhS2E+/ZrDys8wv3VtTXo7uBIWT2Hb9Jm+e+h0WI+AUzSQyLZflg+C/AsUtpxV4VunvW1WwqdAs34f5oVKgSsBkwPa5GQGzEV6XaFj6kFyLggVXzS+NE9plrNnu/9BbIMCHwBsqdh5vbJ5B5a8qw+cKorE3HvzFQCi1N/xGCz4Q2wyA/gDpT+QNEjs/fmu1Pc4tnRs0+Cw/+W+RNpzljlhSrSi34HMB5OpSIPCxkXTjSnzevdwvupESU8UvZg2OFAzoWj1Idms+Mjho69IKTOTXLz54PF4YkxI80HPgcjDeCiEvLELRetNch75+eXglZr9OnYfCFvI1Us18OVH0+xOZkT30Onk0l7Qb8H0Wde8BZrRh5FUC7zBLM3G1AgGyG9K6PCIRn07EHc+XrtemYRUNIXUXVPrjWhVHRIx3B0UVOncbW01mNV+0AWvFDweT4wJqRQIU60FV1x9C2WclXoUZFHFkIWLoug63oSuxyLVAKjms7VWDXrlZhPcYac39t1481bjMfJEOZp6NgDRjEzLvQMXO4wdxNLOaQpcMi0v8mVKwhSoxBWCTmv7vQBUBuKRx3fFjltNUaEIrxQ8Hk+MCa0UOFw1uMzqRctVZeVKuIzapjoMZXBVFMTOAfV/1/Em3Bt44D92Z/lciSqgcnMU5lw1XG1s2KRUA61C9PRfCkvHI2dgE9n59PfsO2ffn680RAVYs2ltvxfXUhm2oDjsXmuxw/xMmkGhatLMi6wBh+9fNNGqiuVM20BB+8wY6URPv3JytQbl2orUaKwCUxEZomxRFtNn+XlItgPmysx54GaBKT2anIqc6PxukMM4Ok934ntUpkajNx88Hk8MrxRS0EfaLAVgGpnTRmvT56swI4B0hWByRqrcCiV7w/ZnK2rjeDQ5C3kwkC3QyBaYVEY90LVQ6vlI4SNFcHWgOxN5hKK+/OhK1Kaudnil4PF4YkxKpVBlW26drDDqrM8A1amCvJhURBeUgqDlyi+g/iHRLauuT81TGH7tWeew5jyYAqB416kL9nq3IXr5fltNBO5UNCkE1zL75BilY3Ud34B3Tv3E6bMuTMpBoZa4Vmri+4zXIKCzbXN77Oeu401hs1xiZPgKhJWXawCPcuw+rEyX+a1rEyW/XZOfygwYJudjFX0pdLjZkNdcIMgB3I0bgeF4JGbHQuDEQHUxC9588Hg8MSa1UiiTI2Ejr3lSdPSulcKgJa2tIGUwPtgyOLljMlECLsWkIKpwQroeK09Kc1o8Qt7uXBTd2NP/7XAbmS1VxCx4peDxeGJMyHoKNqbPmB6LbtTfK4pNdeQ9X5ZPolYKQV+W7Fg4lgiEqbWTsYpCrPw4vNqyLc/CRTVw9eGyf/fQaTzXla4U9PoI+r0mivbxpOXJXbuja334WdWS3qRmXespTBrzgR7AtAeU3nMZGKo4Bj+OqQx9nlDpqmMYaHA49NqesKoTFWWplZOxqsEgjQt9B8JakS5JWFVw8FiyEYvJnLCFMVfR2LdqvPng8XhiTHil4DJrZ+1jUwb6MVz2dT2vK7wydJWoOHxaFryu8uMD7vUS86LXg+w+PIYWTRlwxVCl85FjawdP7/FlSFIGVUQm0rHWrHwjXJ7UvydFYna8UvB4PDFKKQUhxBcA/BFUZ+kDUG3j2gA8CuDXoNrTf0ZKebHkdWaSZwbXyTOj833JL5B27lotieaBR9hFjq+gCvTKKWEUX+v66n0JtbbpgUgNzD/fh+7DKqiHfifT+atQDFRcZd1ioGWJ8hf0QKmCpYgUg+5L4P6DKn0JSomo331np31fFworBSHEIgB/AuBGKeUHAVwB4B4AXwHwVSnlKgCnANxX/jI9Hk+9KOtTaAYwTQhxCcBVAIYA3ALgPwbvPwzgywC+UfI8RmzLj5wqyrYR3Eaz2WtVq4O8pdxMZcNp5orUQe17UNazzyVhq51g2kb728rQpxVypZ4Ra1aq4iykGEyMDe+r2WpD5Ld4pfSxCg8KUspBIcTfATgCVXzux8EVvSulHA12OwZgUemrrIgqHtSshCjXnhT14obN6svaHutIXR/qeb4Nm1aFLeF27lZOt23Lep0+azMpTIMBbettvg6vdqpBYR3ig4OJWi4/UpIULTeXSZAqYz7MAXAngBUAFgKYDuC2HJ+/XwixTwixT46NfwCVx+NRlDEfPgrgLSnlCAAIIR4H8BsAZgshmgO1sBjAoOnDUsodAHYAKqKxyAXkdS7mCTyykXd5sBYKIavOI5kN9TQVxpu2aUoZUORg73BgLlV0fN2UaF/ehJ5+ZS4MBK70rWz5kSIN5+spoDWAZ0wCwMEgS372zA/kXpYssyR5BMAmIcRVQggB4FYABwH8FMAng322A3iqxDk8Hk+dKeNT2COEeAzAqwBGAfwr1Mz/QwCPCiH+Ktj2zSoutArKdJTKKqGWt8BrWfi5uVOR+iVQVeL3g0IA1PIj+QSOBt2lwsrMFWDq8XCh7wDWrFRBXzRTc/JUac6DS+AT7ylJ38UjjscvtfogpfwSgC9pm98EsLHMcWuFydwoMlDoKwGmwaFoVSYTJvNg1VXncVuHEnobb/4wgPfPAJAGOQzDdnTBykA3tsSqKuXB1ASXzjN0fjV6WGMdIB6xaGvu4krRyEdVfxPA4E/qaj54PJ5JyITPfShLGecjVww0m9Oo7Bob4WqekIlw191bATRO85ZGwdTMtqXZvQCKjt7GDgA23hyPdyDHJpBUA1UsP+ZRCWSqRA5H9bsfzNcnCIBXCh6PR+N9qxSqjDjkWYx8xie1kDcHgo7F1YFXBvkhP4JyQGa3rrcxfUEHhs6vBhAVonn+qWdKHbMIFKQExJ2baU5Nrjbr4mh8P5HVfbpsFd25c+aFg8C169XK+h//l98pdUxPhC302RTJSGYDr5pED+TAD36g3myuJt3cZibwQYAwrXTo28qUf/fmg8fjiTEplUIe04BHfKU1leXvZZHX0UjH3ba53SuDEuRN0zY1g+GQQuCQCfd8EI7X038pUXeRz/plnY0jw28YVYELPMIxr7PRKwWPxxNjUiqFPI69tFToshGJrmnd22+/BYD3H9jQVUCRIC1TRCIFOZn20x13I8Nv4Pmn4p2ZgPxLkaQkGrFgK+GVgsfjiTFplMJ4ljwD8tdW+M21bQC8QnDBpgx4Ydg0v0LLquuBIB+CYyqgoudLUKiyyb635Ta0r3gX3W/NTmw3KYS0Yq7zW9cmMix5jwcbVLq/Y+FY7qzeSTMolCGPqTB3zjzj8mPeWARVSTni+aeeSY1F2P9yX2Jpqqf/krGCcBrccca/zHp8PjUynSiYejyY2sxtvFlFgu79udvgQIwMJ00GomizWJ28kYsuzkeKaOw63sSaAP0/p3N488Hj8cSYNEqhXpWTuUrgJoPpvDYFsrOzGwBwaPseACrtWZ8BDr2m3us7Ny13cFRCMvZ3s3NHm1ddpWbOM/NfB6BmGFIStgawExVSAybHY2/zdeHvbJLppMyoxJ3JFAhbxhtMB07R7EdeubnocmUWXil4PJ4Yk67BLFC8YnNWKLMOD1TKqxT0Y9iuqwhVdJSKwq6V/2PNyikNm4Nhalmvv8/pPhxlUFJFZu5ktM3Ceu0EV1zVgSm02YbtWruON4X+hS9+5X+8vxrMcvKaEvqDmRWv4Fq0wlaNKW/hi7wDlsncyDtQ9J2bpv4PTJ2dnQirF5OEbpRBIm/sAi+60o3swYA/XEVpar0xd56DmzMzPfKxyDV788Hj8cSYlEqB0CMKq3BC8hmeK5JanCvr/OPRR+LF/isBAEeCdX8e/z/RHJLchMiiY+FYIs/BlbxRjPFZX6mHqpY/XfBKwePxxJjUSkHHtfQaz5zUbf+0fIY8eRamQiymYxRpI27D1nqON6IlX4IN7m+Igqg+C6DxFIMpsGn/y31oX67KqZFiyFric1mCtPkMyrSfL+5vyM/7alBwpeqH0ZV6mSCmwcFlIEiDTAo8+a1gS2MMDlnp1FRJKQolTnfYpZkOWfEIWWQ5H8tSxEHqzQePxxPjfacUXNKZTTN03hZ1NkyRkKbj1ytKsyqODKkZiUcDjrdaAMzLlRs2rQpn+eHuRzKP4eokTEtuSsO0H6/MTMrFlK/igl+S9Hg8pclUCkKIbwG4A8BJKeUHg21zAXwPwHIAhwF8Skp5Kugp+TUAtwM4B+APpJSv1ubSq0OfiV1VQdnIybRjFVUl7/7ql7mWKasIcOKEfokgZwMYX6WQFdCkz9Lcn2CKWtQzSMv6E3R4mjYQz4jkiiFvxGNeXJTCPyHZYv4BAC9IKa8B8ELwMwB8HMA1wb/7AXyjmsv0eDz1IlMpSClfEkIs1zbfCeAjweuHAfwMwJ8H2/9ZqoSKl4UQs4UQbVLKoaouuGpcy6ZVQdE8h7SlVH2JczyCmUyQYlh6sgvDr6mViUbpc8mXDmnGpRm6pz+9gEmt6kxkzfq2Qi61ypIs6mhsZQ/6CQC0qLMIwFG237FgW8MOCmUGhCxHYJq5kLdVXda+efMibHCTQjclbO9xopiHalxWvLpSEUxNW3r6L4UPHJfrOrbBgL/nEqeQRtrDneZc1AcS18IrrpRefZBSSiFE7lRLIcT9UCYGmpq8v9PjaRSKDgrDZBYIIdoAnAy2DwJYwvZbHGxLIKXcAWAHoFKnC15H3ckzw9c6CKre5oKrE5KnWtcKW2ASBSXRjHrh6G4MXEz2cSC4E09f7iMFwGsu8lwGUguhalixKvycrhpcuz2Z34tMnegaufljLrxSREUUnaKfBrA9eL0dwFNs++8LxSYApxvZn+DxeJK4LEl+F8qpOE8IcQzAlwD8DYDvCyHuAzAA4FPB7s9CLUf2QS1J/mENrnnCkuW/0LtLuXSWqgVlC7RQLkT7rVsrcTCaCrHmoWXJFkCbLc+e7MKFqfSTUhHzW9cmgpRott8/nDzu2PC+cLstuMlW/NWVaGnyQKJcXvx1Uo3kVWwuqw+fTnnrVsO+EsDncl3BBMO1HVzaZ9IoWnRFJ+9gUUWFplrjOghEqdDZ6/jTF3Rg4KJ6vWyqGsRGhqegVatqTQ/9yPAbTpGEpjRpl9WNNGiA5dW4aYDg5d/pfPMRXS+dO29cg/fweTyeGO+73IeiFIk4zLvcWXR5tIwpYUunzgstRU5fcJN1v7JLjCZMNRfJuWiSz8umdrGajOr/+a1rE0uQFJHZ/ZY9m5HeC2dox+s+e7LL2MyWoPdsuQ8jw2+kns9V4XC8UvB4PDG8UnDE5PyzOQRts77r/rbjV+1otPWVyKsibI4tV/+Abb8LfZFt3tt8XbRdUwiEsZvW4uRx26b1AiiWq6E7E3v6DyRqMGzdcr21n4TO9AUdRqei7fy0X5kCLF4peDyeGBNaKVSRWZh1bH22zlp9qEIhuGAq4DreRV1N5F0+zKsQSB0cPBbtN32B+j9Pr01A+SXmn+9LfT9v8JGuHtasnJJQDyPDSVVlC6YyMb91bahA9Ia0RTIqJ/SgwHMIqh4g8jzcVTsUbYVXTJiWJ009J6rMkTBBiVCP/PwyAOBexPsr6LgMGHwgoFZvLUu2hA/O2ZPRvjaHHbHOYDZw9IdIT2fmZMUd6O+roilx88HsQExWcDbFOvBiLAQNDkUrTwPefPB4PBqTpm1crUuoVUFRNVNVOTZbezwXFVFk2fKWle8BiDpK2ZQDEC0trh5VDW97m69LdSCmYTIbSCHwZUpTu7g0eP6EaYbOez0mVWNyKuoKxXRO27Lm1i3Xh8e44/6/d2ob55WCx+OJMWmUAlFL5+N4UYVS4D4Fvq0KXBQEBTbd1tGEjTdvBWDu0ESzN4c7EYF034E+I1Mh2aVtTcbPUHgz0bJkizEwSL8urhj0mdsUllwlWVmPpDaefFx18Lrr7q3he65KYUI7Gk1MloGganj9Rt7opujAUDQCcuBiB1qCwcA0ABD6QMBJW1WgQcD1MweD/7nzUZfr3HSJHnj1nunhpMGAS3puFujxBCPDyV4T/Fr1waWIczNvQpQ3HzweT4xJoxTyLuNVec5andc1ejIvVZZtyyrLtrRNzTvrFtP804ULQcE+mxrg2BQAZ++gukcbF8Xv0ZGhsfA6TITmgOXYppnfNmvzGd6c4hz9nCd7sut4E2aMdAJA+Dvxc0WqoHihG68UPB5PjAmtFMbbqVgmujEvpoCpvKqhFsFL75x6O6EWyKm4tK0pESxUK3Vg2qYrhrLYZl2bg9FU28AUFEXHoN/9zPzN6ED8PpBK4PthaA/bI64QiuRATKhBYTxMhDIUvd5at4njg0OtoxxdBwHCZTDoOzfNmsBlO67JjKCH8dWTUTwFwa/fpf6hrSgKb2DLTYbQHAkqQNFDzlu+HXqNP/jp0H5rVm7N2DMdbz54PJ4YE0Yp8NmzERuu2trI51UMLmnaVdwDHrdQRjGkOR35rG9z9JnUAd9fl+R9nd25r1E/l+l6jgyNAZ1qCZIUgy3egKsEvhSpkzcCkqphA8qxCAAznD5purb8CVFeKXg8nhgNpRTKzH6N5m8wOQLzOkZr/bvUqlBLHyLFQE5Hm6+ARxza4v/5TOuibIo4HOk6p4dZmPYIRV0FmPYztZM3vW+COxZtUIZq1J3L7fgmvFLweDwxGkoplKFRFAKnyoCjosccj2IrfGVg7yn1vynIadvmdgD2YiKmMOC87B08G56fKxfdr2DyM9h8CrY8BGMo9spkNmZVPSBNCoEgv4QrDTEoXNHUFHZ/zvsANeJgkEU9r9m1n0StBw/ujLx62QYAwJqVVwBwX0vP++U2nZ9MG/4Q8QefHmaXxCZe+9HmwOT7u0QY0u/ZsXAsdDrq53FlfutazBh5KNdnvPng8XhiuLSN+xaAOwCclFJ+MNj2twC2AbgIoB/AH0op3w3eexDAfQAuA/gTKeWPss5xeWwsoRJcZHItZ9y08463MqnSFKmHaUGynZsUUVDOFbmORZ/bOZDv2ufOmReeP7yOqxy7e1kUAzcRSCFkBUmFAUqG49Bnw+XHhdHSZHj+oT04M38zAGBd857Y54Bkc98iNRpdlMI/AbhN2/Y8gA9KKdsB9AJ4EACEEOsA3APguuAz/1MIke8v7/F4xhWXXpIvCSGWa9t+zH58GcAng9d3AnhUSvkegLeEEH0ANgJwW1dhjJdCyIJ8H3mohcOxCqqq/qw7EbkqyBuOrOcGVMGqq84bHJ7nwxn22gXRvrYgJBO6ja/7ABL7DMXDlU2KgrZ1HW+KhTrTe0eG1ON0MFAMM9AZqgdouRJFHJlVOBo/C+B7wetFUIMEcSzYVpp6DwS2BzlvTITtGPr58hyjSlwckmkDhx7RaDIZCJOX3NTaLEv25qsped6a4m0zETJbumkPOT+W68BC51iK+P4zRjpxaES9psHm4OhNmBHMsTyGIXy9MN6yzzVnglNqUBBC/CWAUQDfKfDZ+wHcDwBNTd7f6fE0CoUHBSHEH0A5IG+VUaHHQQBL2G6Lg20JpJQ7AOwAVI1GwBwFWA+FYGv0UlQx2Gb2Ms1qx9MEyTI33FrPRUqBR/oVcYjlgRQKRf4BZqmvk7U0qTsYw7Tn0SasK/h0mZYdwxk/NBNqR6EpWghxG4A/A/DbUspz7K2nAdwjhLhSCLECwDUA9pa/TI/HUy9cliS/C+AjAOYJIY4B+BLUasOVAJ4XQgDAy1LK/ySlfF0I8X2ompijAD4npbzsejFFnHi1Yjz6MuShSKBXlRTNquw7Nw1nWHAOYA/q4Y4yU/CSTb2Ylh9NPg3Xyssm1cCXCjkzRjoxPVAi5CswKYAjQ2Ph+2nH0o+rOydNxy0TKemy+vBpw+ZvWvb/awB/XeRiGmVA4INTmQevaFfqtGsq8rla47qCwZ2RJwb2AwAOjagH9Nr1N1UW7psH/SF3MSPSsJkiWVGIdP68TsGDo+qcHevHEuct4mAkvIfP4/HEaIjch0bBVhilFklNZT+T18FZa0zLmqmSHlEewlJL2zMORfCdwJXW8+vn3Dt4Nkyf1h2DQP7u1C5LmHmPmQdSCLQM2YXN6Fhovn9FStd5peDxeGK875VC2vJgLZVBWo6Hi4KoorRbPbG1qqMZ7MX+ebhFd7YFpEYGOpyTnytcihyKHI0ux7NlPdoau/LtWVGOeVTFGbYkSb/TtoXJ34PeM1XbzsIrBY/HE2PSNZjNop4zZpWZnlWEVjcKus0/d868WK8IHX1G3zt4tlRYNuBeos2mFDhFG8raVIKp8K1J3fBCr7TqYOqH0X3g9fdng1kTVT8grpI87UGuxwPbaE5I/oDq8v6dU29HdR2H0isIRdGIbgOjbbnU5Hzk+QqupkqYKp1i/nBMlZqmL+hIrfbME7XoukxVrs+e7AqvI+89MuHNB4/HE6MhzAchxAjU0JZv7aQ2zIO/Do6/jjgT+TqWSSnnZ+3UEIMCAAgh9rnYO/46/HX466jtdXjzwePxxPCDgsfjidFIg8KO8b6AAH8dcfx1xJn019EwPgWPx9MYNJJS8Hg8DUBDDApCiNuEED1CiD4hxAN1OucSIcRPhRAHhRCvCyE+H2yfK4R4Xgjxb8H/c+p0PVcIIf5VCPFM8PMKIcSe4J58TwgxtQ7XMFsI8ZgQ4g0hxCEhxObxuB9CiC8Ef5NfCCG+K4Roqdf9EEJ8SwhxUgjxC7bNeA+E4r8H19QthLihxtfxt8HfplsI8YQQYjZ778HgOnqEEL9V5tzjPigEfSG+DuDjANYB+HTQP6LWjAL4opRyHYBNAD4XnPcBAC9IKa8B8ELwcz34PIBD7OevAPiqlHIVgFNQDXZqzdcAPCelXAtgfXA9db0fQohFAP4EwI1B86EroHqJ1Ot+/BOSfU7S7sHHoUoOXgNVhPgbNb6O+vRbkVKO6z8AmwH8iP38IIAHx+E6ngLwMQA9ANqCbW0Aeupw7sVQX7ZbADwDQEAFpjSb7lGNrmEWgLcQ+JnY9rreD6iWAEcBzIUKw38GwG/V834AWA7gF1n3AMD/AvBp0361uA7tvd8B8J3gdeyZAfAjAJuLnnfclQKiLwFRWa8IV4JmNx8CsAdAq5RyKHjrBIDWOlzCP0AVwqWA+18D8K6UcjT4uR73ZAWAEQDfDsyYh4QQ01Hn+yGlHATwdwCOABgCcBrAK6j//eCk3YPx/O5+FsC/1OI6GmFQGFeEEDMA/ADAn0opY6l3Ug27NV2eEUJQn85XankeB5oB3ADgG1LKD0GFncdMhTrdjzlQncZWAFgIYDqSMnrcqMc9yKJMvxUXGmFQcO4VUTVCiClQA8J3pJSPB5uHhRBtwfttAE7W+DJ+A8BvCyEOA3gUyoT4GoDZQgjKYq3HPTkG4JiUkip+PgY1SNT7fnwUwFtSyhEp5SUAj0Pdo3rfD07aPaj7d5f1W/m9YICq/DoaYVDoAnBN4F2eCuUwebrWJxWqNv03ARySUv49e+tpANuD19uhfA01Q0r5oJRysZRyOdTv/qKU8vcA/BRRj856XMcJAEeFEGuCTbdCleqv6/2AMhs2CSGuCv5GdB11vR8aaffgaQC/H6xCbAJwmpkZlVO3fiu1dBrlcKjcDuVN7Qfwl3U651YoGdgNYH/w73Yoe/4FAP8G4CcA5tbxPnwEwDPB618P/rB9AP4vgCvrcP4NAPYF9+RJAHPG434A+K8A3gDwCwD/G6rHSF3uB4DvQvkyLkGpp/vS7gGUQ/jrwff2ANSKSS2vow/Kd0Df139k+/9lcB09AD5e5tw+otHj8cRoBPPB4/E0EH5Q8Hg8Mfyg4PF4YvhBwePxxPCDgsfjieEHBY/HE8MPCh6PJ4YfFDweT4z/D4BVFV0ilP3rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xad2696ac>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# map each pixel to the centroid value\n",
    "X_recovered = centroids[idx.astype(int),:]\n",
    "X_recovered.shape\n",
    "\n",
    "# reshape to the original dimensions\n",
    "X_recovered = np.reshape(X_recovered, (A.shape[0], A.shape[1], A.shape[2]))\n",
    "X_recovered.shape\n",
    "plt.imshow(X_recovered)\n",
    "X_recovered.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Group C\n",
    "## Principal component analysis\n",
    "\n",
    "PCA is a linear transformation that finds the \"principal components\", or directions of greatest variance, in a data \n",
    "set. It can be used for dimension reduction among other things. In this exercise we're first tasked with implementing \n",
    "PCA and applying it to a simple 2-dimensional data set to see how it works. Let's start off by loading and visualizing\n",
    "the data set.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "from scipy.io import loadmat\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load the data \n",
    "data = loadmat('data/GroupC.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = data['X']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0xacfbbfcc>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArwAAAHVCAYAAAATqShMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAG5NJREFUeJzt3V+Inel9H/Dv09GYnGzdDK3V4JmNs3s1Nxb1GOFSFExr40xCjBGiFw6kF+nF3pTgEJiwKoWSiyLBQEmvAovd4BLXIXXkufC2nhjk4BoaB63H7iS256LGxj5qujJlSGwOjTp9eqEZydrMn3Ok8+99zucDYrXvntX57c6B832f9/d7nlJrDQAAtOpvzboAAACYJIEXAICmCbwAADRN4AUAoGkCLwAATRN4AQBomsALAEDTBF4AAJom8AIA0LRLk/hD3/GOd9SXXnppEn80AAAkSd54440f1FovX/S6iQTel156Kffu3ZvEHw0AAEmSUsp3h3mdlgYAAJom8AIA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaQIvAABNE3gBAGjapVkXAABwlp29frZ3D3L/cJDVlV62NtdzfWNt1mXRMQIvADCXdvb6uXlnP4OHR0mS/uEgN+/sJ4nQy0i0NAAAc2l79+Bx2D0xeHiU7d2DGVVEVwm8AMBcun84GOk6nEXgBQDm0upKb6TrcBaBFwCYS1ub6+ktLz11rbe8lK3N9RlVRFcZWgMA5tLJYJpdGnheAi8AMLeub6wJuDw3LQ0AADRN4AUAoGkCLwAATRN4AQBomsALAEDTBF4AAJom8AIA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaQIvAABNE3gBAGiawAsAQNMEXgAAmibwAgDQNIEXAICmXRh4SynrpZSv/divvyyl/Po0igMAgOd16aIX1FoPkrwnSUopS0n6ST474boAAGAsRm1p+GCS/1Fr/e4kigEAgHEbNfB+NMmnT/sHpZRXSin3Sin3Hjx48PyVAQDAGAwdeEspb0vykST/6bR/Xmt9rdZ6tdZ69fLly+OqDwAAnssoK7y/mOSrtdb/NaliAABg3EYJvL+cM9oZAABgXg0VeEspLyT5UJI7ky0HAADG68JtyZKk1vqjJH9vwrUAAMDYOWkNAICmCbwAADRN4AUAoGkCLwAATRN4AQBomsALAEDTBF4AAJom8AIA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaZdmXQAAAOOzs9fP9u5B7h8OsrrSy9bmeq5vrHX2fcZB4AUAaMTOXj837+xn8PAoSdI/HOTmnf0kGWsYndb7jIuWBgCARmzvHjwOoScGD4+yvXvQyfcZF4EXAKAR9w8HI12f9/cZF4EXAKARqyu9ka7P+/uMi8ALANCIrc319JaXnrrWW17K1uZ6J99nXAytAQA04mRgbNK7J0zrfcal1FrH/odevXq13rt3b+x/LgAAnCilvFFrvXrR67Q0AADQNIEXAICmCbwAADRN4AUAoGkCLwAATRN4AQBomsALAEDTBF4AAJom8AIA0DRHCwMwkp29fmeOEwVIBF4ARrCz18/NO/sZPDxKkvQPB7l5Zz9JhF5gbmlpAGBo27sHj8PuicHDo2zvHsyoIoCLCbwADO3+4WCk6wDzQOAFYGirK72RrgPMA4EXgKFtba6nt7z01LXe8lK2NtdnVBHAxQytATC0k8E0uzQAXSLwAjCS6xtrAi7QKVoaAABomsALAEDTBF4AAJom8AIA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE1z8AQAQEfs7PWddPgMBF4AgA7Y2evn5p39DB4eJUn6h4PcvLOfJELvBbQ0AAB0wPbuweOwe2Lw8Cjbuwczqqg7BF4AgA64fzgY6TpPCLwAAB2wutIb6TpPCLwAsGB29vq5dvtuXn719Vy7fTc7e/1Zl8QQtjbX01teeupab3kpW5vrM6qoOwytAcACMfjUXSc/H7s0jE7gBYAFct7gk+A0/65vrPk5PQMtDQCwQAw+sYgEXgBYIAafWEQCLwAsEINPLCI9vACwQAw+sYgEXgBYMAafWDRaGgAAaJrACwBA0wReAACaJvACANA0Q2sAAB2xs9e3w8YzEHgBADpgZ6+fm3f2Hx8N3T8c5Oad/SQRei+gpQEAoAO2dw8eh90Tg4dH2d49mFFF3WGFFwCYax7jP3L/cDDSdZ6wwgsAzK2Tx/j9w0FqnjzG39nrz7q0qVtd6Y10nScEXgBgbnmM/8TW5np6y0tPXestL2Vrc31GFXWHlgYAYG55jP/ESRuH9o7RCbwAwNxaXemlf0q4XdTH+Nc31gTcZ6ClAQCew85eP9du383Lr76ea7fvLmRv6SR5jM84WOEFgGdkX9TJ8xifcRgq8JZSVpJ8PMm7k9Qk/7zW+t8mWRgAzLvzBqoEsvHxGJ/nNewK779L8vla6z8tpbwtyU9OsCYA6AQDVdANF/bwllJ+Ksn7k3wiSWqtf11rPZx0YQAw7+yLCt0wzNDay0keJPndUspeKeXjpZQX3vqiUsorpZR7pZR7Dx48GHuhADBvDFRBNwwTeC8leW+S36m1biT5UZJX3/qiWutrtdartdarly9fHnOZADB/rm+s5daNK1lb6aUkWVvp5daNKwvXb2qnCubdMD2830/y/VrrV47//jM5JfACwCJa9IEqO1XQBReu8NZa/yLJ90opJ89nPpjkGxOtCgDoBEf/0gXD7tLwa0k+dbxDw7eT/OrkSgIAusJOFXTBUIG31vq1JFcnXAsA0DGO/qULHC0MADwzO1XQBY4WBjjFzl7fUaYwBEf/0gUCL8BbmDpnnnTh5mvRd6pg/mlpAHgLU+fMi5Obr/7hIDVPbr7scwujscIL8BamzpmGYVZuz7v5sqIKw7PCC/AWZ02XmzpnXIZduXXzBeMh8AK8halzJm3Ythk3XzAeAi/AW1zfWMutG1eyttJLSbK20sutG1c8QmZshl25dfMF46GHF+AUps6ZpGEPa7DlF4yHwAsAU7a1uf7U1nfJ2Su3br7g+Qm8ADBlVm5hugReAJgBK7cwPYbWAABomsALAEDTBF4AAJom8AIA0DSBFwCApgm8AAA0TeAFAKBp9uEFAJiwnb2+g0ZmSOAFAJignb3+U0dJ9w8HuXlnP0mE3ikReAGgAVYQ59f27sHjsHti8PAo27sHfkZTIvACQMdZQZxv9w8HI11n/AytAUDHnbeCyOytrvRGus74CbwAMAd29vq5dvtuXn719Vy7fTc7e/2h/10riPNta3M9veWlp671lpeytbk+o4oWj5YGAJix521JWF3ppX9KuLWCOB9OfoZ6rGdH4AWAGXveoaatzfWnAnNiBXHeXN9YE3BnSOAFgBl73pYEK4hwPoEXAGZsHC0JVhDhbIbWAGDGDDXBZFnhBYAZ05IAkyXwAsAc0JIAk6OlAQCAplnhBQDgQjt7/c623Qi8AACc63kPR5k1LQ0AAJzrvMNRukDgBQDgXM97OMqsCbwAAJzrrENQRjkcZZYEXgC4wM5eP9du383Lr76ea7fvZmevP+uSYKq6fjiKoTUAOEfXh3VgHLp+OIrACwDnOG9Ypytf9jAOXT4cRUsDAJyj68M6gMALAOfq+rAOIPACMCGtDHp1fVgH0MMLwAS0NOg1zWGdLh/dCvNM4AVg7Fob9JrGsE5LNwkwb7Q0ADB2Br1G1/WjW2GeCbwAjJ1Br9G5SYDJEXgBGDuDXqNzkwCTI/ACMHbXN9Zy68aVrK30UpKsrfRy68YVvajncJMAk2NoDYCJ6PKpTLPQ9aNbYZ4JvAAwJ9wkwGRoaQAAoGkCLwAATRN4AQBomsALAEDTBF4AAJom8AIA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaZdmXQAAzIudvX62dw9y/3CQ1ZVetjbXc31jbdZlwcQsymde4AVgaGd9Obbwpbmz18/NO/sZPDxKkvQPB7l5Zz9JOvffAsNYpM+8lgYAhnLy5dg/HKTmyZfjv9rZP/X6zl5/1iWPZHv34PEX/4nBw6Ns7x7MqCKYrEX6zAu8AAzlrC/HT3/le018ad4/HIx0HbpukT7zWhp4Li08xgSGc9aX4FGtI71+Xq2u9NI/pebVld4MqoHJW6TPvBVentlZjze79hgTGM5ZX4JLpYz0+nm1tbme3vLSU9d6y0vZ2lyfUUUwWYv0mRd4eWaL1PsDnP3l+Mv/8Gea+NK8vrGWWzeuZG2ll5JkbaWXWzeueGpFsxbpM6+lgWe2SL0/wJOp7dPamK7+7N9tor3p+sZaJ+uGZ7Uon3mBl2e2SL0/wCNnfTkuypcm0E1aGnhmi9T7AwB011ArvKWU7yT5qyRHSf5vrfXqJIuiG857vAkAMC9GaWn4J7XWH0ysEjrJY0wAYN5paQAAoGnDBt6a5I9KKW+UUl457QWllFdKKfdKKfcePHgwvgoBAOA5DNvS8HO11n4p5e8n+UIp5Vu11i/9+Atqra8leS1Jrl69evqxOwAwIU5+BM4yVOCttfaP//pmKeWzSd6X5Evn/1sAsyUALY6Tkx9PDsM5OfkxiZ85cHFLQynlhVLK209+n+Tnk/zZpAsDeB6Ovl4sTn4EzjNMD+9PJ/lyKeXrSf40yeu11s9PtiyA5yMALRYnPwLnubClodb67ST/YAq1AIyNALRYnPwInMe2ZECTzgo6AlCbnPwInEfgBZokAC2W6xtruXXjStZWeilJ1lZ6uXXjioE1IMloJ60BdIajrxePkx+Bswi8QLMEIAASLQ0AADRO4AUAoGkCLwAATdPDCzzmKF4AWiTwAkmeHMV7cjrZyVG8SYTeKXPjATBeWhqAJI7inRcnNx79w0Fqntx47Oz1Z10aQGcJvEASR/HOCzceAOMn8AJJHMU7L9x4AIyfwAskcRTvvHDjATB+Ai+Q5NFg2q0bV7K20ktJsrbSy60bVwxLTZkbD4Dxs0sD8JijeGfv5P+/XRoAxkfgBZgzbjwAxktLAwAATRN4AQBompYGABaKk+xg8Qi8ACwMR2gzDDdF7dHSAMDCcJIdF3G8d5us8AKwMJxkt5hGWbE976bIKm93CbwALIzVlV76p4RbJ9l131mhdtQ2FjdFbdLSAMDCcJJdm85rQxi1jcXx3m0SeIEz7ez1c+323bz86uu5dvuuHjY6zxHabTov1I66YuumqE1aGoBTmWanVU6ya895oXbUNhbHe7dJ4AVOZXAD6IrzQu3W5vpTN+/JxSu2borao6UBOJXBDaArzmtD0MZCYoUXOINpdqArLmpDsGKLwAuc6lkeA7bMyUsw34RaziPwAqcyuPGEAT6AbhN4gTNZMXnEAB9AtxlaA7iAAT6AbhN4AS7g5CWAbhN4AS7g5KXROKEPmDd6eAEuYIBveAb8gHkk8AIMwQDfcAz4AfNISwMAY2PAD5hHAi8AY2PAD5hHAi8AY2PAD5hHeniBznPs7/ww4AfMI4EX6DS7AswfA37AvNHSAHTaebsCAEBihRc6y2P8R+wKAMBFrPBCB508xu8fDlLz5DH+Ip5oZVcAAC4i8EIHeYz/hF0BALiIlgboII/xn7ArAAAXEXihg1ZXeumfEm4X9TG+XQEAOI+WBuggj/EBYHhWeKGDPMYHgOEJvNBRHuMDwHC0NAAA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaZdmXQAwv3b2+tnePcj9w0FWV3rZ2lzP9Y21WZcFACMReIFT7ez1c/POfgYPj5Ik/cNBbt7ZTxKhF4BO0dIAnGp79+Bx2D0xeHiU7d2DGVUEAM9G4AVOdf9wMNJ1AJhXAi9wqtWV3kjXAWBeCbzAqbY219NbXnrqWm95KVub6zOqCACejaE14FQng2l2aQCg6wRe4EzXN9YEXAA6T0sDAABNE3gBAGiawAsAQNMEXgAAmibwAgDQNIEXAICmCbwAADRN4AUAoGlDB95SylIpZa+U8rlJFgQAAOM0ygrvx5J8c1KFAADAJAx1tHAp5cUkv5Tk3yT5jYlWBCyMnb1+tncPcv9wkNWVXrY21x1lDMDYDRV4k/x2kt9M8vazXlBKeSXJK0nyrne96/krA5q2s9fPzTv7GTw8SpL0Dwe5eWc/SYReAMbqwpaGUsqHk7xZa33jvNfVWl+rtV6ttV69fPny2AoE2rS9e/A47J4YPDzK9u7BjCoCoFXD9PBeS/KRUsp3kvx+kg+UUn5volUBzbt/OBjpOgA8qwsDb631Zq31xVrrS0k+muRurfVXJl4Z0LTVld5I1wHgWdmHF5iJrc319JaXnrrWW17K1ub6jCoCoFXDDq0lSWqtf5zkjydSCbBQTgbT7NIAwKSNFHgBxun6xpqAC8DECbwwQ/ahBYDJE3hhRuxDCwDTIfDCjJy3D63A+4RVcACel8ALM2If2otZBQdgHGxLBjNiH9qLOY0NgHEQeGFG7EN7MavgAIyDwAszcn1jLbduXMnaSi8lydpKL7duXPGo/sdYBQdgHPTwwgzZh/Z8W5vrT/XwJlbBARidwAvMLaexATAOAi8w16yCA/C8BF7gb7D3LQAtEXiBp9j7FoDW2KUBeIq9bwFojcALPMXetwC0RksDNGCcPberK730Twm39r4FoKus8ELHnfTc9g8HqXnSc7uz1z/z9ddu383Lr76ea7fv/o3XOQEOgNYIvNBxo/TcDhOOnQAHQGu0NEDHjdJze144/vFAa+9bAFpihRc67qze2tOuG0gDYBEJvNBxo/TcjhKOAaAVAi903Cg9twbSAFhEenihAcP23J68xrHBACwSgRcWjIG06Rnn/sgAPDuBF2ACTraAO9kV42QLuCRCL8CU6eEFmIBR9kcGYLIEXoAJsAUcwPwQeAEmwBZwAPND4AWYAFvAAcwPQ2sAE2ALOID5IfACTIgt4ADmg5YGAACaJvACANA0gRcAgKYJvAAANE3gBQCgaQIvAABNE3gBAGiawAsAQNMEXgAAmibwAgDQNIEXAICmCbwAADRN4AUAoGkCLwAATbs06wIAumJnr5/t3YPcPxxkdaWXrc31XN9Ym3VZAFxA4AUYws5ePzfv7Gfw8ChJ0j8c5Oad/SQRegHmnJYGgCFs7x48DrsnBg+Psr17MKOKABiWwAswhPuHg5GuAzA/BF6AIayu9Ea6DsD8EHgBhrC1uZ7e8tJT13rLS9naXJ9RRQAMy9AawBBOBtPs0gDQPQIvMBEtbuF1fWOt8/8NAItI4AXGzhZeAMwTPbzA2NnCC4B5IvACY2cLLwDmicALjJ0tvACYJwIv2dnr59rtu3n51ddz7fbd7Oz1Z10SHWcLLwDmiaG1BWe4iEmwhRcA80TgXXDnDRcJJzwPW3gBMC+0NCw4w0UAQOsE3gVnuAgAaJ3Au+AMFwEArdPDu+AMFwEArRN4MVwEADRNSwMAAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaQIvAABNE3gBAGiawAsAQNMEXgAAmibwAgDQNIEXAICmXRh4Syk/UUr501LK10spf15K+a1pFAYAAONwaYjX/J8kH6i1/rCUspzky6WU/1Jr/ZMJ1wYAAM/twsBba61Jfnj8t8vHv+okiwIAgHEZqoe3lLJUSvlakjeTfKHW+pVTXvNKKeVeKeXegwcPxl0nAAA8k6ECb631qNb6niQvJnlfKeXdp7zmtVrr1Vrr1cuXL4+7TgAAeCYj7dJQaz1M8sUkvzCZcgAAYLyG2aXhcill5fj3vSQfSvKtSRcGAADjMMwuDe9M8slSylIeBeQ/qLV+brJlAQDAeAyzS8N/T7IxhVoAAGDsnLQGAEDThmlp6ISdvX62dw9y/3CQ1ZVetjbXc31jbdZlAQAwY00E3p29fm7e2c/g4VGSpH84yM07+0ki9AIALLgmWhq2dw8eh90Tg4dH2d49mFFFAADMiyYC7/3DwUjXAQBYHE0E3tWV3kjXAQBYHE0E3q3N9fSWl5661lteytbm+owqAgBgXjQxtHYymGaXBgAA3qqJwJs8Cr0CLgAAb9VESwMAAJxF4AUAoGkCLwAATRN4AQBomsALAEDTBF4AAJom8AIA0DSBFwCApgm8AAA0TeAFAKBpAi8AAE0TeAEAaJrACwBA0wReAACaJvACANA0gRcAgKaVWuv4/9BSHiT57hn/+B1JfjD2N6VLfAYWm58/PgOLzc+fcX4GfrbWevmiF00k8J77hqXcq7VeneqbMld8Bhabnz8+A4vNz59ZfAa0NAAA0DSBFwCAps0i8L42g/dkvvgMLDY/f3wGFpufP1P/DEy9hxcAAKZJSwMAAE0TeAEAaNrUAm8p5d+XUt4spfzZtN6T+VFK+ZlSyhdLKd8opfx5KeVjs66J6Sql/EQp5U9LKV8//gz81qxrYvpKKUullL1SyudmXQvTV0r5Tillv5TytVLKvVnXw3SVUlZKKZ8ppXyrlPLNUso/mtp7T6uHt5Ty/iQ/TPIfaq3vnsqbMjdKKe9M8s5a61dLKW9P8kaS67XWb8y4NKaklFKSvFBr/WEpZTnJl5N8rNb6JzMujSkqpfxGkqtJ/k6t9cOzrofpKqV8J8nVWquDJxZQKeWTSf5rrfXjpZS3JfnJWuvhNN57aiu8tdYvJfnf03o/5kut9X/WWr96/Pu/SvLNJGuzrYppqo/88Phvl49/mZpdIKWUF5P8UpKPz7oWYLpKKT+V5P1JPpEktda/nlbYTfTwMgOllJeSbCT5ymwrYdqOH2d/LcmbSb5Qa/UZWCy/neQ3k/y/WRfCzNQkf1RKeaOU8sqsi2GqXk7yIMnvHrc1fbyU8sK03lzgZapKKX87yR8m+fVa61/Ouh6mq9Z6VGt9T5IXk7yvlKK9aUGUUj6c5M1a6xuzroWZ+rla63uT/GKSf3Hc7shiuJTkvUl+p9a6keRHSV6d1psLvEzNcd/mHyb5VK31zqzrYXaOH2N9MckvzLoWpuZako8c93D+fpIPlFJ+b7YlMW211v7xX99M8tkk75ttRUzR95N8/8ee7H0mjwLwVAi8TMXxwNInknyz1vpvZ10P01dKuVxKWTn+fS/Jh5J8a7ZVMS211pu11hdrrS8l+WiSu7XWX5lxWUxRKeWF46HlHD/K/vkkdm5aELXWv0jyvVLK+vGlDyaZ2uD6pWm9USnl00n+cZJ3lFK+n+Rf11o/Ma33Z+auJflnSfaPeziT5F/WWv/zDGtiut6Z5JOllKU8utn+g1qrralgcfx0ks8+Wv/IpST/sdb6+dmWxJT9WpJPHe/Q8O0kvzqtN3a0MAAATdPSAABA0wReAACaJvACANA0gRcAgKYJvAAANE3gBQCgaQIvAABN+//hoNdcnjbtJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xace6256c>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter(X[:, 0], X[:, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Here we difine the function for PCA. Write a few lines of code to \n",
    "* normalize the data\n",
    "* Computing the covariance matrix\n",
    "* Performing SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2)"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.array([[2,4],[1,3],[0,0],[0,0]])\n",
    "U, S, V = np.linalg.svd(a)\n",
    "V.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2,)"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pca(X):\n",
    "    # normalize the features\n",
    "    X = (X-np.mean(X,axis=0))\n",
    "    \n",
    "    # compute the covariance matrix\n",
    "    X = np.matrix(X)\n",
    "    cov = np.cov(X,rowvar=False)\n",
    "    \n",
    "    # perform SVD\n",
    "    U, S, V = np.linalg.svd(cov)\n",
    "    \n",
    "    return U, S, V\n",
    "\n",
    "U, S, V = pca(X)\n",
    "\n",
    "c=np.dot(np.dot(U,S),V)\n",
    "c.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-14-ea74567147f4>, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-14-ea74567147f4>\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m    Now that we have the principal components (matrix U), we can use these to project the original data into a lower-\u001b[0m\n\u001b[0m           ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "Now that we have the principal components (matrix U), we can use these to project the original data into a lower-\n",
    "dimensional space. For this task write a line in the a function that computes the projection and **selects only the top K \n",
    "components**, effectively reducing the number of dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 2)"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def project_data(X, U, k):\n",
    "    U_reduced = U[:k,:]\n",
    "    return np.dot(X, U_reduced.T)\n",
    "Z = project_data(X, U, 1)\n",
    "Z\n",
    "def recover_data(X, U, k):\n",
    "    U_reduced = U[:,:k]\n",
    "    return np.dot(X, U_reduced.T)\n",
    "X_recovered = recover_data(Z, U, 1)\n",
    "U_reduced.T.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0xa796ab8c>"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAHVCAYAAADywj0dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3W9spelZHvDrjtcthxBw20zpjnfDRi1yJdgmg6wkaBCCRWD+RMtoG9QgoA2CLiDahpYa7SBEk6jVLjWitI1EtE1UBRIgNEzcJQFMpI1UiJpFnvUuJtlYCjQlOZM2Q4Lzp3HprHn6YexhxmuPj2fs855j/36S5eP3PD7n/nC0e83r+7mfaq0FAABOuud1XQAAAIwCwRgAACIYAwBAEsEYAACSCMYAAJBEMAYAgCSCMQAAJBGMAQAgiWAMAABJkju6euMXvvCF7Z577unq7QEAOCEuXrz4p621U/ut6ywY33PPPVleXu7q7QEAOCGq6n8Osk4rBQAARDAGAIAkgjEAACQRjAEAIIlgDAAASQRjAABIIhgDAEASwRgAAJIIxgAAkEQwBgCAJIIxAAAkEYwBACCJYAwAAEkEYwAASCIYAwBAkuSOrgsAAOD4WlzpZ2FpLZfWN3J6qpf5uZmcOzPddVm7EowBADgSiyv9nL+wmo0rm0mS/vpGzl9YTZKRDMdaKQAAOBILS2vXQvG2jSubWVha66iimxOMAQA4EpfWNw50vWuCMQAAR+L0VO9A17smGAMAcCTm52bSm5y44VpvciLzczMdVXRzNt8BAHAktjfYmUoBAMCJd+7M9MgG4Z20UgAAQARjAABIIhgDAEASPcYAANzEOB3pfLsEYwAAnmNxpZ/XPfbBrG9cuXZt1I90vl1aKQAAuMHiSj/nL6zeEIq3jfKRzrdroDvGVfXRJJ9Lspnk2dba7I7nvyHJf03yP7YuXWitveHwygQAYFgWltaycWVzz+dH9Ujn23WQVopvbK396U2e/93W2itvtyAAALq1X/Ad1SOdb5dWCgAAbnCz4DvKRzrfrkGDcUvyO1V1saoe3GPN11bV01X1W1X1VbstqKoHq2q5qpYvX758SwUDAHC05udm0puceM71v/bFk3n4gXuP5ca7ZPBWiq9rrfWr6m8meW9Vfbi19t+ue/7JJF/RWvt8VX17ksUkX7nzRVprjyZ5NElmZ2fbbdYOAMAR2A6+J2VM27Zq7WD5tKpel+TzrbWfvcmajyaZvVlP8uzsbFteXj7QewMAwEFV1cWdwyN2s28rRVU9v6pesP04ybck+cMda/5WVdXW45dtve6nbqVwAADowiCtFF+e5F1bufeOJL/cWvvtqvrhJGmtvSnJq5L8SFU9m2QjyavbQW9FAwBAh/YNxq21P07ykl2uv+m6x29M8sbDLQ0AAIbHuDYAAIhgDAAASQRjAABIIhgDAEASwRgAAJIIxgAAkGTwI6EBAOjY4kr/xB3TPEyCMQDAGFhc6ef8hdVsXNlMkvTXN3L+wmqSCMeHRCsFAMAYWFhauxaKt21c2czC0lpHFR0/gjEAwBi4tL5xoOscnGAMADAGTk/1DnSdgxOMAQDGwPzcTHqTEzdc601OZH5upqOKjh+b7wAAxsD2BjtTKY6OYAwAMCbOnZkWhI+QYAwA0BFziUeLYAwA0AFziUePzXcAAB0wl3j0CMYAAB0wl3j0CMYAAB0wl3j0CMYAAB0wl3j02HwHANABc4lHj2AMANARc4lHi1YKAACIYAwAAEkEYwAASCIYAwBAEsEYAACSmEoBADCwxZW+8WrHmGAMADCAxZV+zl9YzcaVzSRJf30j5y+sJolwfExopQAAGMDC0tq1ULxt48pmFpbWOqqIwyYYAwAM4NL6xoGuM34EYwCAAZye6h3oOuNHMAYAGMD83Ex6kxM3XOtNTmR+bqajijhsNt8BAAxge4OdqRTHl2AMADCgc2emBeFjTCsFAABEMAYAgCSCMQAAJBGMAQAgic13AMAxt7jSN0mCgQjGAMCxtbjSz/kLq9eOcu6vb+T8hdUkEY55Dq0UAMCxtbC0di0Ub9u4spmFpbWOKmKUDXTHuKo+muRzSTaTPNtam93xfCX590m+PckXkrymtfbk4ZYKADCY7faJ/vrGrs9f2uM6J9tBWim+sbX2p3s8921JvnLr6+VJfmHrOwDAUO1sn9jN6aneECtiXBxWK8V3JvnFdtUHkkxV1Z2H9NoAAAPbrX3ier3JiczPzQyxIsbFoMG4JfmdqrpYVQ/u8vx0ko9d9/PHt67doKoerKrlqlq+fPnywasFANjHzdokpqd6efiBe228Y1eDtlJ8XWutX1V/M8l7q+rDrbX/dtA3a609muTRJJmdnW0H/X0AgP2cnurt2ls8PdXL+x+6r4OKGBcD3TFurfW3vn8yybuSvGzHkn6Su6/7+a6tawAAQzU/N5Pe5MQN17RPMIh9g3FVPb+qXrD9OMm3JPnDHcseS/IP66pXJPlMa+0Th14tAMA+zp2ZzsMP3JvpqV4q2icY3CCtFF+e5F1XJ7LljiS/3Fr77ar64SRprb0pyW/m6qi2j+TquLbvP5pyAYCT6iAn2J07My0Ic2D7BuPW2h8necku19903eOW5EcPtzQAgKucYMcwOPkOABh5TrBjGA5ywAcAwFDsbJtwgh3DIBgDACNlt7aJytVDFXZygh2HSSsFADBSdmubaElqxzoj2DhsgjEAMFL2ao9oiRFsHCmtFADASHFyHV1xxxgA6MTiSj9nH3k8L37oPTn7yONZXLl6aK6T6+iKO8YAwNANMpd40MM84LAIxgDA0N1sLvH2qXWCMMOmlQIAGLq9NtiZS0yXBGMAYOj2mj9sLjFdEowBgKGzwY5RpMcYABg6G+wYRYIxANAJG+wYNVopAAAggjEAACQRjAEAIIlgDAAASQRjAABIIhgDAEASwRgAAJIIxgAAkEQwBgCAJIIxAAAkEYwBACBJckfXBQAA3Vpc6WdhaS2X1jdyeqqX+bmZnDsz3XVZMHSCMQCcYIsr/Zy/sJqNK5tJkv76Rs5fWE0S4ZgTRysFAJxgC0tr10Lxto0rm1lYWuuoIuiOO8YAcILsbJvor2/suu7SHtfhOBOMAeCE2K1topK0XdaenuoNtTYYBYIxABxzP7W4ml954mPZbM+NwC15TjjuTU5kfm5mWOXByNBjDADH2E8truZtH/iTXUPxtpZkeqqX2vr+8AP32njHieSOMQAcY7/yxMf2XTM91cv7H7pvCNXAaHPHGACOsZvdKU60TcD13DEGgGNsomrPcDztMA+4gWAMAGPuZifXfffL787bPvAnz/md733Fi/Kvz9077FJhpAnGADDG9ju5bjv8bk+lmKjKd7/8bqEYdlFtn96jozI7O9uWl5c7eW8AOC7OPvL4rod02FAHf6mqLrbWZvdb544xAIyRxZV+Xv8bH8yffeHKTdc5uQ4OTjAGgDGxuNLP/DufzpXN/f/a6+Q6ODjj2gBgTCwsrQ0Uio1gg1szcDCuqomqWqmqd+/y3Guq6nJVPbX19YOHWyYAsF97hJPr4PYcpJXitUmeSfKlezz/jtbaP7n9kgCA3Zye6u260S6x2Q4Ow0B3jKvqriTfkeTNR1sOAJxsiyv9nH3k8bz4offk7COPZ3Glf+25+bmZTE7Uc35n8nmldQIOwaCtFD+f5CeS/MVN1vz9qvqDqnpnVd19+6UBwMmyPZO4v76Rlr+cSbwdjs+dmc7Cq16Sv/bFk9d+Z6o3mYXveonWCTgE+7ZSVNUrk3yytXaxqr5hj2W/keRXWmt/XlU/lOStSZ7z95yqejDJg0nyohe96JaLBoDjaGFp7dpBHds2rmxmYWntWvA9d2ZaCIYjMkiP8dkk91fVtyf5oiRfWlVva6197/aC1tqnrlv/5iT/drcXaq09muTR5OoBH7dcNQAcE9cf57zX/xjNJIbh2DcYt9bOJzmfJFt3jP/l9aF46/qdrbVPbP14f65u0gMAbuKnFlfz9g/8yZ6BeJuZxDAct3zAR1W9Iclya+2xJP+squ5P8mySTyd5zeGUBwDH0+JKf6BQbCYxDE+11k1Hw+zsbFteXu7kvQGga2cfeXzP0WvJ1ZnEp6d6mZ+b0VMMt6mqLrbWZvdb50hoAOjAzfqGzSSGbjgSGgA6sFffcCVaJ6AjgjEAdGB+bia9yYkbrlWS73nFi7ROQEe0UgBAB7bD7/aoNv3E0D3BGAAO0fVzifcLuw7rgNEiGAPAIdk+0nn79LrtI52TCMAwBvQYA8AhudmRzsDoE4wB4JDsNYLNkc4wHgRjADgke41gc6QzjAfBGAAOyW4j2BzpDOPD5jsAOCRGsMF4E4wB4BAZwQbjSysFAABEMAYAgCSCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgjGAACQxMl3AJwQiyt9RzUDNyUYA3DsLa70c/7CajaubCZJ+usbOX9hNUmEY+AarRQAHHsLS2vXQvG2jSubWVha66giYBQJxgAce5fWNw50HTiZBGMAjr3TU70DXQdOJsEYgGNhcaWfs488nhc/9J6cfeTxLK70rz03PzeT3uTEDet7kxOZn5sZdpnACLP5DoCxt9/muu0NdqZSADcjGAMw9m62uW47/F4fkAF2o5UCgLFncx1wGARjAMaezXXAYRCMARh7NtcBh0GPMQBjz+Y64DAIxgAcCzbXAbdLMAZgJC2u9N0BBoZKMAZg5Ow3lxjgKNh8B8DIudlcYoCj4o4xAJ3b2TbRN5cY6IBgDECndmubqCRtl7XmEgNHSSsFAJ3arW2iJakd68wlBo6aYAxAp/Zqj2hJpqd6qa3vDz9wr413wJHSSgHAUCyu9PP63/hg/uwLV5IkU73JvO7+r9qzp3h6qpf3P3TfsMsETrCB7xhX1URVrVTVu3d57q9W1Tuq6iNV9URV3XOYRQIw3hZX+pl/59PXQnGSrG9cyfx/eTrf+HdPOc4ZGAkHaaV4bZJn9njuB5L8WWvt7yT5d0l+5nYLA+B4WFzp58d/7elc2Xzudrorf9Hyvg9fzsMP3KttAujcQK0UVXVXku9I8m+S/ItdlnxnktdtPX5nkjdWVbXWdttUDMAJsT1xYvMm/zu4tL7hOGdgJAx6x/jnk/xEkr/Y4/npJB9Lktbas0k+k+Rv7FxUVQ9W1XJVLV++fPkWygVgnOw2cWInI9iAUbFvMK6qVyb5ZGvt4u2+WWvt0dbabGtt9tSpU7f7cgCMuP0O5Jh8XuklBkbGIHeMzya5v6o+muRXk9xXVW/bsaaf5O4kqao7knxZkk8dYp0AjKGb3Q2e6k1m4bteooUCGBn7BuPW2vnW2l2ttXuSvDrJ4621792x7LEk/2jr8au21ugvBjjh5udmdp048fP/4KV56l99i1AMjJRbnmNcVW9IstxaeyzJW5L8UlV9JMmnczVAA3DCbQffhaW1XFrfyOmpXubnZgRiYCRVVzd2Z2dn2/LycifvDQDAyVFVF1trs/utcyQ0AABEMAYAgCSCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgjGAACQRDAGAIAkgjEAACQRjAEAIIlgDAAASZI7ui4AgKOzuNLPwtJaLq1v5PRUL/NzMzl3ZrrrsgBGkmAMcEwtrvRz/sJqNq5sJkn66xs5f2E1SYRjgF1opQA4phaW1q6F4m0bVzazsLTWUUUAo00wBjimLq1vHOg6wEknGAMcU6enege6DnDSCcYAx9T83Ex6kxM3XOtNTmR+bqajigBGm813AMfU9gY7UykABiMYA4yhQcewnTszLQgDDEgwBhgzxrABHA09xgBjxhg2gKMhGAOMGWPYAI6GYAwwZoxhAzgagjHAmDGGDeBo2HwHMGaMYQM4GoIxwBgyhg3g8GmlAACACMYAAJBEMAYAgCSCMQAAJLH5DmCoFlf6pkkAjCjBGGBIFlf6OX9h9dpxzv31jZy/sJokwjHACNBKATAkC0tr10Lxto0rm1lYWuuoIgCuJxgDDMml9Y0DXQdguARjgCE5PdU70HUAhkswBhiS+bmZ9CYnbrjWm5zI/NxMRxUBcD2b7wCGZHuDnakUAKNJMAYYonNnpgVhgBElGAPcIjOJAY4XwRjgFphJDHD87Lv5rqq+qKp+v6qerqoPVtXrd1nzmqq6XFVPbX394NGUCzAazCQGOH4GuWP850nua619vqomk/xeVf1Wa+0DO9a9o7X2Tw6/RIDRYyYxwPGz7x3jdtXnt36c3PpqR1oVwIgzkxjg+BlojnFVTVTVU0k+meS9rbUndln296vqD6rqnVV19x6v82BVLVfV8uXLl2+jbIBumUkMcPwMFIxba5uttZcmuSvJy6rqq3cs+Y0k97TW/l6S9yZ56x6v82hrbba1Nnvq1KnbqRugU+fOTOfhB+7N9FQvlWR6qpeHH7jXxjuAMVatHawroqp+OskXWms/u8fzE0k+3Vr7spu9zuzsbFteXj7QewMAwEFV1cXW2ux+6waZSnGqqqa2HveSfHOSD+9Yc+d1P96f5JmDlQsAAN0aZCrFnUneunUn+HlJfq219u6qekOS5dbaY0n+WVXdn+TZJJ9O8pqjKhgAAI7CgVspDotWCgAAhuHQWikAAOAkEIwBACCCMQAAJBGMAQAgyWBTKQDG1uJKPwtLa7m0vpHTU73Mz804hAOAXQnGwLG1uNLP+Qur2biymSTpr2/k/IXVJBGOAXgOwRg4Vq6/Q/y8qmzuGEm5cWUzC0trgjEAzyEYA8fGzjvEO0PxtkvrG8MsC4AxYfMdcGwsLK1dC8U3c3qqN4RqABg3gjFwbAxyJ7g3OZH5uZkhVAPAuBGMgWNjrzvBE1WpJNNTvTz8wL36iwHYlR5j4NiYn5u5occ4uXqHWBgGYBCCMTA29ptJvP3Y3GIAboVgDIyFQWcSnzszLQgDcEv0GANjYbeJE9sziQHgMAjGwFjYa+KEmcQAHBbBGBgLe02cMJMYgMMiGANjYX5uJr3JiRuumUkMwGGy+Q4YCyZOAHDUBGNgbJg4AcBR0koBAAARjAEAIIlWCmCI9ju5DgC6JBgDQzHoyXUA0BWtFMBQOLkOgFEnGAND4eQ6AEadYAwMhZPrABh1gjEwFE6uA2DU2XwHDIWT6wAYdYIxMDROrgNglAnGwC0zlxiA40QwBm7J9/yn/573/9Gnr/1sLjEA487mO+DAdobibeYSAzDOBGPgQBZX+ruG4m3mEgMwrgRj4ED2uyNsLjEA40owBg5kvzvC5hIDMK4EY+BAbnZH+Ozf/us23gEwtgRj4AaLK/2cfeTxvPih9+TsI49ncaV/w/O7nWCXXA3Fb//HXzusMgHg0BnXBlyzuNLP+Qur2biymWT3EWxOsAPguBKMgWsWltauheJt2yPYrg++TrAD4DgSjOEE23lyXX+PjXVGsAFwEgjGcELt1jZRSdoua41gA+Ak2HfzXVV9UVX9flU9XVUfrKrX77Lmr1bVO6rqI1X1RFXdcxTFAodnt7aJlqR2rOtNThjBBsCJMMhUij9Pcl9r7SVJXprkW6vqFTvW/ECSP2ut/Z0k/y7JzxxumcBh26s9oiWZnuqltr4//MC9+okBOBH2baVorbUkn9/6cXLra+dfW78zyeu2Hr8zyRurqrZ+FxhBe/UUT0/18v6H7uugIgDo1kBzjKtqoqqeSvLJJO9trT2xY8l0ko8lSWvt2SSfSfI3dnmdB6tquaqWL1++fHuVA7dlt3nE2iYAOMkGCsattc3W2kuT3JXkZVX11bfyZq21R1trs6212VOnTt3KSwCH5NyZ6Tz8wL3aJgBgy4GmUrTW1qvqfUm+NckfXvdUP8ndST5eVXck+bIknzq0KoEjYR4xAPylQaZSnKqqqa3HvSTfnOTDO5Y9luQfbT1+VZLH9RcDADBOBrljfGeSt1bVRK4G6V9rrb27qt6QZLm19liStyT5par6SJJPJ3n1kVUM7GnngR2OagaAwQ0yleIPkpzZ5fpPX/f4/yb5rsMtDTiI3Q7sOH9hNUmEYwAYwECb74DRt9uBHRtXNrOwtNZRRQAwXgRjOCb2OrBjr+sAwI0EYzgmTk/1DnQdALiRYAxjYnGln7OPPJ4XP/SenH3k8Syu9G943oEdAHB7DjTHGOjGIBvrtr+bSgEAt0YwhjFws4111wdfB3YAwK3TSgFjwMY6ADh67hjDiNntkI7TU730dwnBNtYBwOFxxxhGyE8truafv+Op9Nc30vKXvcTf+HdP2VgHAEdMMIYRsbjSz9s/8CdpO65vXNnM+z58OQ8/cG+mp3qpJNNTvTz8wL36iQHgEGmlgBGxsLT2nFC87dL6ho11AHDEBGPo0PX9xHuF4kQvMQAMg2AMHdk5m3gvleglBoAh0GMMHdltNvFOleR7XvEiLRQAMATuGENHbjaDuBIn1wHAkAnG0JG9ZhNPT/Xy/ofu66AiADjZtFJAR+bnZswmBoAR4o4xdGS7RWLnKXdaJwCgG4IxdMhsYgAYHVopAAAggjEAACQRjAEAIIlgDAAASQRjAABIIhgDAEASwRgAAJIIxgAAkEQwBgCAJIIxAAAkEYwBACCJYAwAAEkEYwAASCIYAwBAEsEYAACSCMYAAJBEMAYAgCSCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgwQjKvq7qp6X1V9qKo+WFWv3WXNN1TVZ6rqqa2vnz6acgEA4GjcMcCaZ5P8eGvtyap6QZKLVfXe1tqHdqz73dbaKw+/RAAAOHr73jFurX2itfbk1uPPJXkmyfRRFwYAAMN0oB7jqronyZkkT+zy9NdW1dNV9VtV9VV7/P6DVbVcVcuXL18+cLEAAHBUBg7GVfUlSX49yY+11j674+knk3xFa+0lSf5jksXdXqO19mhrbba1Nnvq1KlbrRkAAA7dQMG4qiZzNRS/vbV2YefzrbXPttY+v/X4N5NMVtULD7VSAAA4QoNMpagkb0nyTGvt5/ZY87e21qWqXrb1up86zEIBAOAoDTKV4myS70uyWlVPbV37ySQvSpLW2puSvCrJj1TVs0k2kry6tdaOoF4AADgS+wbj1trvJal91rwxyRsPqygAABg2J98BAEAEYwAASCIYAwBAEsEYAACSCMYAAJBEMAYAgCSCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgjGAACQRDAGAIAkgjEAACQRjAEAIIlgDAAASQRjAABIIhgDAEASwRgAAJIkd3RdAKNhcaWfhaW1XFrfyOmpXubnZnLuzHTXZQEADI1gTBZX+jl/YTUbVzaTJP31jZy/sJokwjEAcGJopSALS2vXQvG2jSubWVha66giAIDhE4zJpfWNA10HADiOBGNyeqp3oOsAAMeRYEzm52bSm5y44VpvciLzczMdVQQAMHw233Ftg52pFADASSYYk+RqOBaEAYCTTCsFAABEMAYAgCSCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgjGAACQRDAGAIAkgjEAACQRjAEAIIlgDAAASQYIxlV1d1W9r6o+VFUfrKrX7rKmquo/VNVHquoPquprjqZcAAA4GncMsObZJD/eWnuyql6Q5GJVvbe19qHr1nxbkq/c+np5kl/Y+g4AAGNh3zvGrbVPtNae3Hr8uSTPJJnesew7k/xiu+oDSaaq6s5DrxYAAI7IgXqMq+qeJGeSPLHjqekkH7vu54/nueEZAABG1sDBuKq+JMmvJ/mx1tpnb+XNqurBqlququXLly/fyksAAMCRGCgYV9Vkrobit7fWLuyypJ/k7ut+vmvr2g1aa4+21mZba7OnTp26lXoBAOBIDDKVopK8JckzrbWf22PZY0n+4dZ0ilck+Uxr7ROHWCcAABypQaZSnE3yfUlWq+qprWs/meRFSdJae1OS30zy7Uk+kuQLSb7/8EsFAICjs28wbq39XpLaZ01L8qOHVRQAAAybk+8AACCCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgjGAACQRDAGAIAkgjEAACQRjAEAIIlgDAAASQRjAABIIhgDAEASwRgAAJIIxgAAkCS5o+sChmlxpZ+FpbVcWt/I6ale5udmcu7MdNdlAQAwAk5MMF5c6ef8hdVsXNlMkvTXN3L+wmqSCMcAAJycVoqFpbVroXjbxpXNLCytdVQRAACj5MQE40vrGwe6DgDAyXJigvHpqd6BrgMAcLKcmGA8PzeT3uTEDdd6kxOZn5vpqCIAAEbJidl8t73BzlQKAAB2c2KCcXI1HAvCAADs5sS0UgAAwM0IxgAAEMEYAACSCMYAAJBEMAYAgCSCMQAAJBGMAQAgiWAMAABJBGMAAEgiGAMAQBLBGAAAkgjGAACQRDAGAIAkgjEAACQRjAEAIIlgDAAASZJqrXXzxlWXk/zP23iJFyb500Mqh/Hms0Dic8BVPgckPgdcdf3n4Ctaa6f2+4XOgvHtqqrl1tps13XQPZ8FEp8DrvI5IPE54Kpb+RxopQAAgAjGAACQZLyD8aNdF8DI8Fkg8TngKp8DEp8Drjrw52Bse4wBAOAwjfMdYwAAODSCMQAAZAyDcVXdXVXvq6oPVdUHq+q1XdfE8FXVF1XV71fV01ufg9d3XRPdqaqJqlqpqnd3XQvdqaqPVtVqVT1VVctd10M3qmqqqt5ZVR+uqmeq6mu7ronhqqqZrf8ObH99tqp+bKDfHbce46q6M8mdrbUnq+oFSS4mOdda+1DHpTFEVVVJnt9a+3xVTSb5vSSvba19oOPS6EBV/Ysks0m+tLX2yq7roRtV9dEks601BzucYFX11iS/21p7c1X9lSRf3Fpb77ouulFVE0n6SV7eWtv3YLmxu2PcWvtEa+3JrcefS/JMkuluq2LY2lWf3/pxcutrvP6Vx6GoqruSfEeSN3ddC9CtqvqyJF+f5C1J0lr7f0LxifdNSf5okFCcjGEwvl5V3ZPkTJInuq2ELmz9+fypJJ9M8t7Wms/ByfTzSX4iyV90XQida0l+p6ouVtWDXRdDJ16c5HKS/7zVXvXmqnp+10XRqVcn+ZVBF49tMK6qL0ny60l+rLV53SfoAAABmElEQVT22a7rYfhaa5uttZcmuSvJy6rqq7uuieGqqlcm+WRr7WLXtTASvq619jVJvi3Jj1bV13ddEEN3R5KvSfILrbUzSf5Pkoe6LYmubLXS3J/kvwz6O2MZjLd6Sn89ydtbaxe6rodubf2Z7H1JvrXrWhi6s0nu3+ot/dUk91XV27otia601vpb3z+Z5F1JXtZtRXTg40k+ft1fEN+Zq0GZk+nbkjzZWvvfg/7C2AXjrU1Xb0nyTGvt57quh25U1amqmtp63EvyzUk+3G1VDFtr7Xxr7a7W2j25+ueyx1tr39txWXSgqp6/tSE7W386/5Ykf9htVQxba+1/JflYVc1sXfqmJDbnn1zfnQO0USRX/+Qwbs4m+b4kq1v9pUnyk6213+ywJobvziRv3dpt+rwkv9ZaM6oLTq4vT/Kuq/dOckeSX26t/Xa3JdGRf5rk7Vt/Rv/jJN/fcT10YOsfyN+c5IcO9HvjNq4NAACOwti1UgAAwFEQjAEAIIIxAAAkEYwBACCJYAwAAEkEYwAASCIYAwBAkuT/A/SyKDzz+EvxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa7bb10ec>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,8))\n",
    "ax.scatter([X_recovered[:, 0]],[X_recovered[:, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Group D\n",
    "## Application of Principal component analysis\n",
    "\n",
    "PCA is a linear transformation that finds the \"principal components\", or directions of greatest variance, in a data \n",
    "set. It can be used for dimension reduction among other things. In this exercise we're first tasked with implementing \n",
    "PCA and applying it to a simple 2-dimensional data set to see how it works.  Our task here is to apply PCA to images of faces.  By using the same dimension reduction \n",
    "techniques we can capture the \"essence\" of the images using much less data than the original images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "from scipy.io import loadmat\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the face image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'hdf5storage'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-179-64dcf4e825c7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0mhdf5storage\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mfaces\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhdf5storage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mloadmat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data/GroupD.mat'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m#faces = loadmat('data/GroupD.mat')\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mX\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfaces\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'X'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'hdf5storage'"
     ]
    }
   ],
   "source": [
    "import hdf5storage\n",
    "faces = hdf5storage.loadmat('data/GroupD.mat')\n",
    "#faces = loadmat('data/GroupD.mat')\n",
    "X = faces['X']\n",
    "X.shape\n",
    "face = np.reshape(X[3,:], (32, 32))\n",
    "plt.imshow(face)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Here we difine the function for PCA. Write a few lines of code to \n",
    "* normalize the data\n",
    "* Computing the covariance matrix\n",
    "* Performing SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-19-856ac270d390>, line 3)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-19-856ac270d390>\"\u001b[0;36m, line \u001b[0;32m3\u001b[0m\n\u001b[0;31m    X =\u001b[0m\n\u001b[0m        ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "def pca(X):\n",
    "    # normalize the features\n",
    "    X = \n",
    "    \n",
    "    # compute the covariance matrix\n",
    "    X = np.matrix(X)\n",
    "    cov = \n",
    "    \n",
    "    # perform SVD\n",
    "    U, S, V = \n",
    "    \n",
    "    return U, S, V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'pca' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-20-f831ecaabb3f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mV\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpca\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mV\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'pca' is not defined"
     ]
    }
   ],
   "source": [
    "U, S, V = pca(X)\n",
    "U, S, V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the principal components (matrix U), we can use these to project the original data into a lower-\n",
    "dimensional space. For this task write a line in the a function that computes the projection and **selects only the top K \n",
    "components**, effectively reducing the number of dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-21-976583ad5767>, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-21-976583ad5767>\"\u001b[0;36m, line \u001b[0;32m2\u001b[0m\n\u001b[0;31m    U_reduced =\u001b[0m\n\u001b[0m                ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "def project_data(X, U, k):\n",
    "    U_reduced = \n",
    "    return np.dot(X, U_reduced)\n",
    "Z = project_data(X, U, 1)\n",
    "Z\n",
    "def recover_data(Z, U, k):\n",
    "    U_reduced = U[:,:k]\n",
    "    return np.dot(Z, U_reduced.T)\n",
    "X_recovered = recover_data(Z, U, 1)\n",
    "X_recovered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'pca' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-b8495dfb8cbf>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mV\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpca\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mZ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mproject_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'pca' is not defined"
     ]
    }
   ],
   "source": [
    "U, S, V = pca(X)\n",
    "Z = project_data(X, U, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'recover_data' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-0395ad4ac261>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mX_recovered\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrecover_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mZ\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mface\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_recovered\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m32\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m32\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mface\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'recover_data' is not defined"
     ]
    }
   ],
   "source": [
    "X_recovered = recover_data(Z, U, 100)\n",
    "face = np.reshape(X_recovered[3,:], (32, 32))\n",
    "plt.imshow(face)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}