lectures/svd_part2_pca.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis\n",
    "## Intuition\n",
    "Imagine, that you have a dataset of points. Your goal is to choose orthogonal axes, that describe your data the most informative way. To be precise, we choose first axis in such a way, that maximize the variance (expressiveness) of the projected data. All the following axes have to be orthogonal to the previously chosen ones, while satisfy largest possible variance of the projections. \n",
    "\n",
    "Let's take a look at the simple 2d data. We have a set of blue points on the plane. We can easily see that the projections on the first axis (red dots) have maximum variance at the final position of the animation. The second (and the last) axis should  be orthogonal to the previous one.\n",
    "![](https://i.stack.imgur.com/lNHqt.gif)\n",
    "[source](https://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues)\n",
    "\n",
    "This idea could be used in a variety ways. For example, it might happen, that projection of complex data on the principal plane (only 2 components) bring you enough intuition for clustering. The picture below plots projection of the labeled dataset onto the first to principal components (PC's), we can clearly see, that only two vectors (these PC's) would be enogh to differ Finnish people from Italian in particular dataset (celiac disease (Dubois et al. 2010))\n",
    "![](https://nla.skoltech.ru/files/pca_example.png)\n",
    "[source](https://privefl.github.io/bigsnpr/articles/how-to-PCA.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem\n",
    "The first component should be defined in order to maximize variance, i.e. the sum of all squared projections of data points to our vector ${\\mathbf  {w}}_{{(1)}}$, which implies the following optimization problem:\n",
    "$$\n",
    "{\\mathbf  {w}}_{{(1)}}={\\underset  {\\Vert {\\mathbf  {w}}\\Vert =1}{\\operatorname {\\arg \\,max}}}\\,\\left\\{\\sum _{i}\\left({\\mathbf  {a}}_{{(i)}}\\cdot {\\mathbf  {w}}\\right)^{2}\\right\\}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    " \\mathbf {w} _{(1)}={\\underset {\\Vert \\mathbf {w} \\Vert =1}{\\operatorname {\\arg \\,max} }}\\,\\{\\Vert \\mathbf {Aw} \\Vert ^{2}\\}={\\underset {\\Vert \\mathbf {w} \\Vert =1}{\\operatorname {\\arg \\,max} }}\\,\\left\\{\\mathbf {w} ^{\\top}\\mathbf {A^{\\top}} \\mathbf {Aw} \\right\\}\n",
    "$$\n",
    "\n",
    "since we are looking for the unit vector, we can reformulate the problem:\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mathbf {w} _{(1)}={\\operatorname {\\arg \\,max} }\\,\\left\\{{\\frac {\\mathbf {w} ^{\\top}\\mathbf {A^{\\top}} \\mathbf {Aw} }{\\mathbf {w} ^{\\top}\\mathbf {w} }}\\right\\}}\n",
    "$$\n",
    "\n",
    "It is [known](https://en.wikipedia.org/wiki/Rayleigh_quotient), that for positive semidefinite matrix $A^\\top A$ such vector is nothing else, but eigenvector of $A^\\top A$, which corresponds to the largest eigenvalue. The following components will give you the same results (eigenvectors).\n",
    "\n",
    "So, we can conclude, that the following mapping:\n",
    "\n",
    "$$\n",
    "\\underset{n \\times k}{\\Pi} = \\underset{n \\times d}{A} \\cdot \\underset{d \\times k}{W} \n",
    "$$\n",
    "\n",
    "describes the projection of data onto the $k$ principal components, where $W$ contains first (by the size of eigenvalues) $k$ eigenvectors of $A^\\top A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll briefly derive how SVD decomposition could lead us to the PCA.\n",
    "\n",
    "Firstly, we write down SVD decomposition of our matrix:\n",
    "\n",
    "$$\n",
    "A = U \\Sigma W^\\top\n",
    "$$\n",
    "\n",
    "and to its transpose:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "A^\\top\n",
    "&= (U \\Sigma W^\\top)^\\top \\\\\n",
    "&= (W^\\top)^\\top \\Sigma^\\top U^\\top \\\\\n",
    "&= W \\Sigma^\\top U^\\top \\\\\n",
    "&= W \\Sigma U^\\top\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Then, consider matrix $A A^\\top$:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "A^\\top A\n",
    "&= (W \\Sigma U^\\top)(U \\Sigma V^\\top)  \\\\\n",
    "&= W \\Sigma I \\Sigma W^\\top \\\\\n",
    "&= W \\Sigma \\Sigma W^\\top \\\\\n",
    "&= W \\Sigma^2 W^\\top\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Which corresponds to the eigendecomposition of matrix $A^\\top A$, where $W$ stands for the matrix of eigenvectors of $A^\\top A$, while $\\Sigma^2$ contains eigenvalues of $A^\\top A$.\n",
    "\n",
    "At the end:\n",
    "$$\n",
    "\\begin{align}\n",
    "\\Pi &= A \\cdot W =\\\\\n",
    " &= U \\Sigma W^\\top W = U \\Sigma\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The latter formula provide us with easy way to compute PCA via SVD with any number of principal components:\n",
    "$$\n",
    "\\Pi_r = U_r \\Sigma_r\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example. Iris dataset\n",
    "Consider the classical Iris dataset\n",
    "![](https://sebastianraschka.com/images/blog/2015/principal_component_analysis_files/iris.png)\n",
    "[source](https://sebastianraschka.com/Articles/2015_pca_in_3_steps.html)\n",
    "\n",
    "We have the dataset matrix $A \\in \\mathbb{R}^{150 \\times 4}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "🤖: Dataset contains 150 points in 4-dimensional space\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "iris        = load_iris()\n",
    "A           = iris['data']\n",
    "labels      = iris['target']\n",
    "classes     = iris['target_names']\n",
    "label_names = np.array([classes[label] for label in labels])\n",
    "colours     = ['red', 'blue', 'green']\n",
    "print('🤖: Dataset contains {} points in {}-dimensional space'.format(*A.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data normalization with zero mean and unit variance\n",
    "A_std = StandardScaler().fit_transform(A)\n",
    "\n",
    "u,s,v = np.linalg.svd(A_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8leWZ//HPxTKGTRah/BCooRZBIiFAQBgZjVCKK9oWBUcsaAfcELHLWOqrsgz91f7qwmirFJGBTlVAtINQRESJqAWFREVWdYRSFAVElrCVwPX743kSAzwnOeTk5CTh+3698so5z3I/17kJufIs93WbuyMiInKiWqkOQEREqiYlCBERiaQEISIikZQgREQkkhKEiIhEUoIQEZFIShAiIhJJCUJERCIpQYiISKQ6qQ4gEc2bN/f09PTIdfv376dBgwaVG1ANov5LjPovcerDxJTWf3l5eTvdvUVZbVTrBJGens6qVasi1+Xm5pKTk1O5AdUg6r/EqP8Spz5MTGn9Z2Z/i6cNXWISEZFIShAiIhJJCUJERCIpQYiISCQlCBERiZS0BGFm081su5mtKbGsmZm9YmYfhd+bhsvNzB41s4/NbLWZdUtWXCIiEp9knkHMAC47YdnPgVfdvT3wavge4HKgffg1EngiiXGJiEgckpYg3H0ZsOuExdcAM8PXM4FrSyz/owdWAE3MrFWyYhMRkbJV9kC5lu6+DcDdt5nZN8LlrYG/l9hua7hsWyXHJyJVwDNvb2Hm2wd5YuPyVIcSU6ezz2Tc1RmpDiOpqspIaotY5pEbmo0kuAxFy5Ytyc3NjWywoKAg5jopm/ovMeq/xMx8+yBb9h4Fdqc6lJi2HttLbu6OVIcRU0X8DFZ2gvjCzFqFZw+tgO3h8q1A2xLbtQE+i2rA3acCUwGys7M91lByDdNPjPovMeq/xARnDrt5+d7LUx1KtVURP4OV/Zjri8Cw8PUwYF6J5T8Mn2bqBewpuhQlIiKpkbQzCDN7FsgBmpvZVmAc8AAwx8x+BGwBrgs3XwhcAXwMHABuTlZcIiISn6QlCHe/IcaqfhHbOnBnsmIREZFTp5HUIiISSQlCREQiKUGIiEikqjIOQqTGqA6DvKq6ddv2cna9VEchOoMQqWDz3vuULfuOpTqMaq1TqzPpfbb+fk01/QuIJME3G9Vi9q29Ux1GtaaR6KmnMwgREYmkBCEiIpGUIEREJJIShIiIRFKCEBGRSEoQIiISSQlCREQiKUGIiEgkJQgREYmkBCEiIpGUIEREJJIShIiIRFKCEBGRSEoQIiISSQlCREQiKUGIiEgkJQgREYmkBCEiIpGUIEREJJIShIiIRFKCEBGRSEoQIiISqU6qA5Cq55m3tzDz7YM8sXF5qkOpltZt28vZ9VIdhUjidAYhJ5n33qds2Xcs1WFUW51anUnvs/W3l1R/+imWSN9sVIvZt/ZOdRjVVm5ubqpDEEmYziBERCRSShKEmd1jZmvNbI2ZPWtmaWbWzszeNrOPzGy2mf1TKmITEZFApScIM2sNjAay3f0CoDYwBPgN8Ii7twe+An5U2bGJiMjXUnWJqQ5Qz8zqAPWBbUBfYG64fiZwbYpiExERUpAg3P1T4EFgC0Fi2APkAbvdvTDcbCvQurJjExGRr1X6U0xm1hS4BmgH7AaeAy6P2NRj7D8SGAnQsmXLmE+LFBQU6EmSctq9+yBHjx5V/yVAP3+JUx8mpiL6LxWPuX4H2OTuOwDM7AXgn4EmZlYnPItoA3wWtbO7TwWmAmRnZ3tOTk7kQXJzc4m1Tkr3xMbl7N69W/2XAP38JU59mJiK6L9U3IPYAvQys/pmZkA/YB2wFBgUbjMMmJeC2EREJJSKexBvE9yMzgc+CGOYCtwL/NjMPgbOAp6q7NhERORrKRlJ7e7jgHEnLP4E6JmCcEREJIJGUouISCQlCBERiaQEISIikZQgREQkkhKEiIhEUoIQEZFIZSYIM2tpZk+Z2Uvh+05mpkqrIiI1XDxnEDOAl4Gzw/cfAmOSFZCIiFQN8SSI5u4+BzgGENZKOprUqEREJOXiSRD7zewswuqqZtaLoES3iIjUYPGU2vgx8CJwrpm9BbTg66J6IiJSQ5WZINw938wuAToABmx09yNJj0xERFIqnqeY7gQauvtad18DNDSzO5IfmoiIpFI89yBGuPvuojfu/hUwInkhiYhIVRBPgqgVTuwDgJnVBv4peSGJiEhVEM9N6peBOWY2heBJptuARUmNSkREUi6eBHEvcCtwO8FN6sXAtGQGJSIiqRfPU0zHgCfCLxEROU2UmSDM7CJgPHBOuL0B7u7fSm5oIiKSSvFcYnoKuAfIQyU2REROG/EkiD3u/lLSIxERkSolngSx1Mx+C7wAHC5a6O75SYtKRERSLp4EcWH4PbvEMgf6Vnw4IiJSVcTzFNOllRGIiIhULfGcQWBmVwIZQFrRMnefmKygREQk9eIp1jcFGAzcRfCI63UEj7yKiEgNFk8tpn929x8CX7n7BKA30Da5YYmISKrFkyAOht8PmNnZwBGgXfJCEhGRqiCeexALzKwJ8Fsgn+AJJtViEhGp4eJ5iuk/wpfPm9kCIM3dNSe1iEgNFzNBmFlfd3/NzL4fsQ53fyG5oYmISCqVdgZxCfAacHXEOicYWS0iIjVUzATh7uPMrBbwkrvPqcSYRESkCij1KaZwLohRFX1QM2tiZnPNbIOZrTez3mbWzMxeMbOPwu9NK/q4IiISv3gec33FzH5qZm3DX+LNzKxZgsf9T2CRu3cEugDrgZ8Dr7p7e+DV8L2IiKRIPI+53hJ+v7PEMgfKNWGQmZ0JXAwMB3D3fwD/MLNrgJxws5lALsF0pyIikgLxPOZa0YPivgXsAP7LzLoQTER0N9DS3beFx9xmZt+o4OPKaezIkSNs3bqVQ4cOVcrxGjduzPr16yvlWDWV+jAxjRs3ZtOmTbRp04a6deuWq414i/VdAHTi+GJ9fyzXEYNjdgPucve3zew/OYXLSWY2EhgJ0LJlS3JzcyO3KygoiLlOSrd790GOHj1ao/qvYcOGtGzZktatW2NmST/e0aNHqV27dtKPU5OpDxNTWFhIQUEB77//PgUFBeVqI545qccRXPrpBCwELgfeBMqbILYCW9397fD9XIIE8YWZtQrPHloB26N2dvepwFSA7Oxsz8nJiTxIbm4usdZJ6Z7YuJzdu3fXqP5bv349bdq0qZTkALBv3z4aNWpUKceqqdSHidm3bx9t2rShoKCA7OzssneIEM9N6kFAP+Bzd7+Z4KbyGeU6GuDunwN/N7MO4aJ+wDrgRWBYuGwYMK+8xxCJUlnJQaSqSPRnPq5ifeHjroXhDebtlPMGdQl3AU+b2WogC/i/wANAfzP7COgfvhepMT7//HOGDBnCueeeS6dOnbjiiiv48MMPk3rMnJwcVq1aVeo2kydP5sCBA8Xvr7jiCnbv3p3UuE7F8OHDmTt3bqnb3H///SxZsqRCjhdPn5VHRcZYWeK5B7EqLNb3JMEN5QLgnUQO6u7vcfwUpkX6JdKuSFXl7nzve99j2LBhzJo1C4D33nuPL774gvPOOy+lsU2ePJmhQ4dSv359ABYuXJjSeMpj4sSqPX/Z0aNHq3yMUco8g3D3O9x9t7tPIfjLflh4qUlE4rR06VLq1q3LbbfdVrwsKyuLf/mXfyE3N5errrqqePmoUaOYMWMGAOnp6fziF7+gd+/eZGdnk5+fz4ABAzj33HOZMmUKQKn7l3T77beTnZ1NRkYG48aNA+DRRx/ls88+49JLL+XSSy8tPubOnTu59957efzxx4v3Hz9+PA899BAAv/3tb+nRoweZmZnFbZ1o8eLF9O7dm27dunHddddRUFDAnj176NChAxs3bgTghhtu4MknnwSCBwl+8pOf0K1bN/r168fOnTtPanPixIn06NGDCy64gJEjR+LuwPFnGenp6YwbN45u3brRuXNnNmzYAMD+/fu55ZZb6NGjB127dmXevOAq9sGDBxkyZAiZmZkMHjyYgwcPnnTcl156ieuvv774fW5uLldffXXMfi2KY+LEifTp04fnnnvuuBhjfY6cnBzuvfdeevbsyXnnnccbb7wBBAnmpz/9KZ07dyYzM5PHHnsMgLy8PC655BK6d+/OgAED2LZtW+S/RXnFc5N6HjAbmOfumyv06CIpMGH+WtZ9trdC2+x09pmMuzoj5vo1a9bQvXv3crXdtm1bli9fzj333MPw4cN56623OHToEBkZGcclnLL86le/olmzZhw9epR+/fqxevVqRo8ezcMPP8zSpUtp3rz5cdsPGTKEMWPGcMcddwAwZ84cFi1axOLFi/noo4945513cHcGDhzIsmXLuPjii4v33blzJ5MmTWLJkiU0aNCA3/zmNzz88MPcf//9/O53v2P48OHcfffdfPXVV4wYMQIIfoF369aNhx56iIkTJ/LrX/+aqVOnHhfTqFGjuP/++wG46aabWLBgQfEv6pKaN29Ofn4+jz/+OA8++CDTpk3jV7/6FX379mX69Ons3r2bnj178p3vfIc//OEP1K9fn9WrV7N69Wq6det2Unv9+/fn1ltvZf/+/TRo0IDZs2czePDgmP2amZkJQFpaGm+++SYAixYtiutzFBYW8s4777Bw4UImTJjAkiVLmDp1Kps2beLdd9+lTp067Nq1iyNHjnDXXXcxb948WrRowezZs7nvvvuYPn16vD8SZYrnHsTDQB9gnZk9Z2aDzCytrJ1EpGIMHDgQgM6dO3PhhRfSqFEjWrRoQVpa2indK5gzZw7dunWja9eurF27lnXr1pW6fdeuXdm+fTufffYZ77//Pk2bNuWb3/wmixcvZvHixXTt2pVu3bqxYcMGPvroo+P2XbFiBevWreOiiy4iKyuLmTNn8re//Q0Iftl27tyZO++8k2nTvp5aplatWsW/dIcOHcqKFStOimnp0qVceOGFdO7cmddee421a9dGxv797wdFqLt3787mzZuB4IzmgQceICsri5ycHA4dOsSWLVtYtmwZQ4cOBSAzM7P4l3tJderU4bLLLmP+/PkUFhbyl7/8hWuuuabMfi36PKfyOaJiX7JkCbfddht16gR/0zdr1oyNGzeyZs0a+vfvT1ZWFpMmTWLr1q2RxyuveAbKvQ68bma1gb7ACGA6cGaFRiJSSUr7Sz9ZMjIyYt5orVOnDseOHSt+f+JgvjPOCB4arFWrVvHroveFhYVl7g+wadMmHnzwQVauXEnTpk0ZPnx4XIMGBw0axNy5c4tvsENwP2Xs2LHceuutMfdzd/r378+zzz570rpjx46xfv166tWrx65du2jTpk1kGyc+gXPo0CHuuOMOVq1aRdu2bRk/fnzMz1DUT7Vr16awsLA4pueff54OHTqctH08T/sMHjyY3//+9zRr1owePXrQqFGjMvu1QYMGJ7VT1ueIFfuJMbo7GRkZLF++vMzYyyueMwjMrB7wA+A2oAdBKQwRiVPfvn05fPhw8fV2gJUrV/L6669zzjnnsG7dOg4fPsyePXt49dVXT6ntePbfu3cvDRo0oHHjxnzxxRe89NJLxesaNWrEvn37ItseMmQIs2bNYu7cuQwaNAiAAQMGMH369OLBV59++inbtx8/bKlXr1689dZbfPzxxwAcOHCg+ImtRx55hPPPP59nn32WW265hSNHjgBB4ihKos888wy9evU6rs2iX6LNmzenoKCgzCebTjRgwAAee+yx4uv97777LgAXX3wxTz/9NBBcCly9enXk/jk5OeTn5/Pkk08WnxmU1q+xlOdzfPe732XKlCnFCWPXrl106NCBHTt2FCeII0eOxDyjKq947kHMBi4EFgG/B3LDx15FJE5mxp///GfGjBnDAw88QFpaGunp6UyePJm2bdty/fXXk5mZSfv27enatesptR3P/l26dKFr165kZGTwrW99i4suuqh43ciRI7n88stp1aoVS5cuPW6/jIwM9u3bR+vWrWnVqhUQ/LJav349vXv3BoKby3/605/4xje+ro7TokULZsyYwQ033MDhw4cBmDRpEgDTpk3jnXfeoVGjRlx88cVMmjSJCRMm0KBBA9auXUv37t1p3LgxTz311HGxNGnShBEjRtC5c2fS09Pp0aPHKfXTL3/5S8aMGUNmZibuTnp6OgsWLOD222/n5ptvJjMzk6ysLHr27Bm5f+3atbnqqquYMWMGM2fOLLNfYynP5/i3f/s3PvzwQzIzM6lbty4jRoxg1KhRzJ07l9GjR7Nnzx4KCwsZM2YMGRkVd4ZsRdk05gZmlwGvuPvRCjtqBcnOzvZYzytrJHX5Df5DMJL65XsvT3UoFWb9+vWcf/75lXY8jQI+dQ0bNjyuJIT6MDFF/Rf1s29mee5e5vDqeO5BLCprGzk1yXiKpiKt27aXs+ulOgoRSbW4ivXVRI+8ktwRrKV5d8tuduw7nLLjA7RpGjsDdGp1JufXL19xL5HyKm9BOUme0zZBpNIl57VIdQjc07/00bs1qZKriJRPzARhZiePFinB3fMrPhwREakqSjuDeCj8nkZQN+l9wIBM4G2CwXMiIlJDxRwH4e6XuvulwN+Abu6e7e7dga7Ax5UVoIiIpEY8A+U6uvsHRW/cfQ1BiW4ROQUNGzY8pe1LFuF78cUXeeCB0ivgxyonfWIxv1NVVLwvWVTOu+qK5yb1ejObBvwJcGAooIlipVqr6KfYyrrpn6iBAwcW12SKpTqWk45XVf9s1bWcd1niOYO4GVgL3A2MIZj9TeW+RcqpaBDnoEGD6NixIzfeeGNx+YdFixbRsWNH+vTpwwsvvFC8z4wZMxg1ahR79uwhPT29uPbSgQMHaNu2LUeOHDnuL/FY7YwfP54HH3yw+P0FF1xQXBDu2muvpXv37mRkZJxURTVKouW8d+zYcVKbJctgjx49WuW8k1zOuyzxzAdxCJgC/Nzdv+fuj4TLRKSc3n33XSZPnsy6dev45JNPikt4jxgxgvnz5/PGG2/w+eefn7Rf48aN6dKlC6+//joA8+fPZ8CAAdStW7d4m3jaiTJ9+nTy8vJYtWoVjz76KF9++WXMbUuW887Pzyc7O5uHH36Yxo0bF5fznjVrVmQ57/z8fC655BImTJhwUrujRo1i5cqVrFmzhoMHD7JgwYLI4xeV87799tuLE15ROe+VK1eydOlSfvazn7F//36eeOKJ4nLe9913H3l5eSe1179/f1asWMH+/fsBTirnvWrVKlavXs3rr79+XK2monLeRYUM4/kcReW8J0+eXNwHJct5r169mhtvvLG4nPfcuXPJy8vjlltu4b777ov5b5IMZSYIMxsIvEdQiwkzyzKzF5MdmEhN1rNnT9q0aUOtWrXIyspi8+bNbNiwgXbt2tG+fXvMrLgE9YkGDx7M7NmzAZg1a9ZJJaXjbedEjz76KF26dKFXr178/e9/P6mEd0kVUc67aJ6EkkqWwV62bJnKeSe5nHdZ4rkHMQ7oCeRCMF2omaUnLySRmq9k2e6SZZ3jKTs9cOBAxo4dy65du8jLy6Nv374nbROrnVilwXNzc1myZAnLly+nfv36xb9gY6mMct5jx45VOe8kl/MuSzz3IArdfU/SIxE5zXXs2JFNmzbxv//7vwCRv3whuJbfs2dP7r77bq666ipq164ddzvp6enk5wdjXPPz89m0aRMAe/bsoWnTptSvX58NGzZETtZTUkWU8+7T5/ihVCeWwS66hxAvlfOuePGcQawxs38FaptZe2A08NfkhiVy+klLS2Pq1KlceeWVNG/enD59+rBmzZrIbQcPHsx1110XWRKltHZ+8IMf8Mc//pGsrCx69OjBeecFT19ddtllTJkyhczMTDp06HDSXAwnqohy3kWXyYqcWAY7aurP0qicd8WLp9x3feA+4LsEI6lfBv6jKtyoTqTcdyqL9VUF8dRiqknl0lXuO/VOLOddFvVhYiqr3PcBggRRubfPRUQkpeKZUe484KdAesnt3f3kO2MiIjGonHf1E889iOcIxkFMA6rcrHIiIpIc8SSIQnd/IumRiCRZ1KOEIjVZWfeYyxLPY67zzewOM2tlZs2KvhI6qkglS0tL48svv0z4P4xIdeHufPnll6SlpZW7jXjOIIaF339W8tjAt8p9VJFK1qZNG7Zu3RpZ/ycZDh06lNB/TFEfJurQoUM0adIk5kDFeMTzFFO7crcuUkXUrVuXdu0q70c5NzeXrl27VtrxaiL1YWIqov9Km3K0r7u/Zmbfj1rv7i9ELRcRkZqhtDOIS4DXgKsj1jmgBCEiUoPFTBDuPi78rrkfREROQ/HcpMbMrgQygOI7Ru6e0PRJZlYbWAV86u5XmVk7YBbQDMgHbnL3fyRyDBERKb945oOYAgwG7iKoxXQdcE4FHPtujp+69DfAI+7eHvgK+FEFHENERMopnnEQ/+zuPwS+cvcJQG+gbSIHNbM2wJUEo7OxYPRSX6CoLu5M4NpEjiEiIomJJ0EUTeB6wMzOBo4AiT4vOBn4d6Bo5pKzgN3uXhi+3wq0TvAYIiKSgHjuQSwwsybAbwnuDTjhX/7lYWZXAdvdPc/McooWR2waOeTVzEYCIwFatmwZWQ8fgsJgsdYBtD50OO6Ya6Lc3M9KXV9W/0np1H+JUx8mpiL6r8z5II7b2OwMIC2RGebM7NfATUAhwU3vM4E/AwOA/+PuhWbWGxjv7gNKa0vzQZTf6TYfRGVT/yVOfZiY0vov4fkgYg2QC9eVe6Ccu48Fxobt5AA/dfcbzew5YBDBk0zDgFObb1BERCpUaZeYogbIFUnGQLl7gVlmNgl4F3iqgtsXEZFTUNpAuaQPkHP3XCA3fP0JED1ZrIiIVLp4xkGcZWaPmlm+meWZ2X+a2VmVEZyIiKROPI+5zgJ2AD8guEewA5idzKBERCT14nnMtZm7/0eJ95PMTIPYRERquHjOIJaa2RAzqxV+XQ/8JdmBiYhIasWTIG4FngEOh1+zgB+b2T4z25vM4EREJHXimVGuUWUEIiIiVUs8TzH96IT3tc1sXPJCEhGRqiCeS0z9zGyhmbUys87ACkBnFSIiNVw8l5j+1cwGAx8AB4Ab3P2tpEcmIiIpFc8lpvYEk/s8D2wGbjKz+kmOS0REUiyeS0zzgV+6+63AJcBHwMqkRiUiIikXz0C5nu6+F8CD2uAPmdmLyQ1LRERSLeYZhJn9O4C77zWz605YnfRCfiIiklqlXWIaUuL12BPWXZaEWEREpAopLUFYjNdR70VEpIYpLUF4jNdR70VEpIYp7SZ1l7DWkgH1StRdMoK5pOU0pjm9S5/TW6QmKG1GudqVGYiIiFQt8YyDEBGR05AShIiIRFKCEBGRSEoQIiISSQlCREQiKUGIiEgkJQgREYmkBCEiIpGUIEREJJIShIiIRFKCEBGRSEoQIiISSQlCREQiKUGIiEgkJQgREYlU6QnCzNqa2VIzW29ma83s7nB5MzN7xcw+Cr83rezYRETka6k4gygEfuLu5wO9gDvNrBPwc+BVd28PvBq+FxGRFKn0BOHu29w9P3y9D1gPtAauAWaGm80Erq3s2ERE5Gvm7qk7uFk6sAy4ANji7k1KrPvK3U+6zGRmI4GRAC1btuw+a9asyLYLCgpo2LBhzGNv33c4kdCrvW80OqPU9eq/0iXaf1I29WFiSuu/Sy+9NM/ds8tqI+ac1MlmZg2B54Ex7r7XzOLaz92nAlMBsrOzPScnJ3K73NxcYq0DeOSVD08t4Brm+pzzSl2v/itdov0nZVMfJqYi+i8lTzGZWV2C5PC0u78QLv7CzFqF61sB21MRm4iIBFLxFJMBTwHr3f3hEqteBIaFr4cB8yo7NhER+VoqLjFdBNwEfGBm74XLfgE8AMwxsx8BW4DrUhCbiIiEKj1BuPubQKwbDv0qMxYREYlNI6lFRCSSEoSIiERSghARkUhKECIiEkkJQkREIilBiIhIJCUIERGJpAQhIiKRlCBERCSSEoSIiERSghARkUgpmw9C5HR2us+nAXBP/9Ln1JDU0xmEiIhEUoIQEZFIShAiIhJJCUJERCIpQYiISCQlCBERiaQEISIikZQgREQkkhKEiIhEUoIQEZFIShAiIhJJCUJERCIpQYiISCQlCBERiaRy3yJSLZ3uJdMro1y6ziBERCSSEoSIiERSghARkUhKECIiEkkJQkREIlWpBGFml5nZRjP72Mx+nup4REROZ1UmQZhZbeD3wOVAJ+AGM+uU2qhERE5fVSZBAD2Bj939E3f/BzALuCbFMYmInLaqUoJoDfy9xPut4TIREUmBqjSS2iKW+UkbmY0ERoZvC8xsY4z2mgM7Kyi2GufHZW+i/iuF+i9x6sPEJNh/58RzjKqUILYCbUu8bwN8duJG7j4VmFpWY2a2yt2zKy6804v6LzHqv8SpDxNTEf1XlS4xrQTam1k7M/snYAjwYopjEhE5bVWZMwh3LzSzUcDLQG1guruvTXFYIiKnrSqTIADcfSGwsIKaK/MylJRK/ZcY9V/i1IeJSbj/zP2k+8AiIiJV6h6EiIhUITUyQahkR/mZ2XQz225ma1IdS3VkZm3NbKmZrTeztWZ2d6pjqk7MLM3M3jGz98P+m5DqmKojM6ttZu+a2YJE2qlxCUIlOxI2A7gs1UFUY4XAT9z9fKAXcKd+/k7JYaCvu3cBsoDLzKxXimOqju4G1ifaSI1LEKhkR0LcfRmwK9VxVFfuvs3d88PX+wj+k6oiQJw8UBC+rRt+6UbpKTCzNsCVwLRE26qJCUIlO6RKMLN0oCvwdmojqV7CyyPvAduBV9xd/XdqJgP/DhxLtKGamCDiKtkhkkxm1hB4Hhjj7ntTHU914u5H3T2LoJpCTzO7INUxVRdmdhWw3d3zKqK9mpgg4irZIZIsZlaXIDk87e4vpDqe6srddwO56J7YqbgIGGhmmwkur/c1sz+Vt7GamCBUskNSxswMeApY7+4Ppzqe6sbMWphZk/B1PeA7wIbURlV9uPtYd2/j7ukEv/tec/eh5W2vxiUIdy8Eikp2rAfmqGRH/MzsWWA50MHMtprZj1IdUzVzEXATwV9u74VfV6Q6qGqkFbDUzFYT/LH3irsn9KimlJ8BHrbSAAAERUlEQVRGUouISKQadwYhIiIVQwlCREQiKUGIiEgkJQgREYmkBCEiIpGUIKRKMLOj4SOha8zsOTOrH2O7hUXPyZ9i+2eb2dwE4ttsZs3Lu391YWbDzezsVMchVYMShFQVB909y90vAP4B3FZypQVqufsV4QjbU+Lun7n7oIoKtgYbDihBCKAEIVXTG8C3zSw9nFfhcSAfaFv0l3yJdU+G8wYsDkfeYmbfNrMl4ZwC+WZ2brj9mnD9cDObZ2aLwnlDxhUd2Mz+x8zywjZHlhVoOPdIfnisV8NlzcJ2VpvZCjPLDJePN7OZYaybzez7Zvb/zOyDMJa64Xabzew34bwI75jZt8Pl55jZq2G7r5rZN8PlM8zsUTP7q5l9YmaDSsT3MzNbGe4zIVwW2XfhftnA0+HZXD0ze8DM1oX7P1gB/7ZSnbi7vvSV8i+gIPxeB5gH3A6kE1Sk7FViu81A83BdIZAVLp8DDA1fvw18L3ydBtQPt18TLhsObAPOAuoBa4DscF2z8HvR8rNKHveEmFsQVA5ud8K+jwHjwtd9gffC1+OBNwlKWHcBDgCXh+v+DFxb4lj3ha9/CCwIX88HhoWvbwH+J3w9A3iO4A++TgTl7gG+SzAvsYXrFgAXl9F3uSX7AtjI1wNqm6T650RflfulMwipKuqFJZ5XAVsI6hkB/M3dV8TYZ5O7vxe+zgPSzawR0Nrd/wzg7ofc/UDEvq+4+5fufhB4AegTLh9tZu8DKwiKPrYvJeZewDJ33xQeq2gejT7Af4fLXgPOMrPG4bqX3P0I8AFQG1gULv+A4Bd3kWdLfO8dvu4NPBO+/u8SMUOQLI65+zqgZbjsu+HXuwRnYB1LfJ6T+i7i8+0FDgHTzOz7BAlNTiN1Uh2ASOigByWeiwV179hfyj6HS7w+SvBXf1S59ygn1phxM8shKA7X290PmFkuwRlILBbRTtHyWMc7DODux8zsiLsXLT/G8f8fPcbrqDaL2z3h+Ab82t3/cFxwwTwVUX13fOPuhWbWE+hHUPhtFMEZkZwmdAYhNYoHcy9sNbNrAczsjBhPRPUP7xXUA64F3gIaA1+FyaEjwRlCaZYDl5hZu/BYzcLly4Abw2U5wE4/9TkhBpf4vjx8/VeCX9SE7b9ZRhsvA7dYMDcFZtbazL5Rxj77gEbh9g2Bxu6+EBhDMAWonEZ0BiE10U3AH8xsInAEuI6TZ9d6k+AyzbeBZ9x9lZl9ANwWVhLdSHCZKSZ33xHeyH7BzGoRzIDWn+Bew3+F7RwAhpXjM5xhZm8T/BF3Q7hsNDDdzH4G7ABuLiO+xWZ2PrA8PBsrAIYSnDHEMgOYYmYHCeZ1n2dmaQRnI/eU43NINaZqrnLaMbPhBDdiR6U6ligWTPaS7e47Ux2LnN50iUlERCLpDEJERCLpDEJERCIpQYiISCQlCBERiaQEISIikZQgREQkkhKEiIhE+v+5KZdCku0CrgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "total_variance = sum(s)\n",
    "variance_explained = [(i / total_variance)*100 for i in sorted(s, reverse=True)]\n",
    "cumulative_variance_explained = np.cumsum(variance_explained)\n",
    "\n",
    "xs = [0.5 + i for i in range(4)]\n",
    "plt.bar(xs, variance_explained, alpha=0.5, align='center',\n",
    "        label='Individual explained variance')\n",
    "plt.step(xs, cumulative_variance_explained, where='mid',\n",
    "         label='Cumulative explained variance')\n",
    "plt.ylabel('Explained variance')\n",
    "plt.xlabel('Principal components')\n",
    "plt.legend(loc='best')\n",
    "plt.xticks([0,1,2,3,4])\n",
    "plt.grid()\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "rank = 2\n",
    "projections = u[:,:rank] @ np.diag(s[:rank])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "setosa red\n",
      "versicolor blue\n",
      "virginica green\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X94VOWZN/DvnZCYRH4oYJESmdjXWqsExGBF374KgtWl/qhV3nWvERHLZimtlWvrvtbmasF20712dVfdvqJSa6UwbeyLWsvW3XXB5LLu0gpxgVHxB60JBlQkKTFpwEByv3+cmTAzOefMOTNz5pyZ+X6ua66Q+XWeQ2bOfc7z3M/9iKqCiIiozO8GEBFRMDAgEBERAAYEIiKKYUAgIiIADAhERBTDgEBERAB8DAgiUiUiL4vILhF5TUTu8astREQEiF/zEEREAJysqv0iUgHgJQB3qOpvfWkQEVGJG+PXhtWIRP2xXytiN86SIyLyiW8BAQBEpBxAO4CzADykqr8zeU4jgEYAqK6ubjjjjDNGHhseHkZZWeEPg3A/gqMY9gEojv0ohn0AgrEfb7311iFVPS3tE1XV9xuAUwC0Aphh97yGhgZN1NraqsWA+xEcxbAPqsWxH8WwD6rB2A8AO9TBsTgQ4VdVDwNoA3CVz00hIipZfmYZnSYip8T+XQ1gIYA3/GoPEVGp83MMYSqA9bFxhDIAv1DVf/GxPUREJc3PLKPdAGb7tX0iIkoWiDEEIiLyHwMCEREBYEAgIqIYBgQiIgLAgEAlKBIB6uqAsjLjZyTid4uIgsHX0hVE+RaJAI2NwMCA8Xtnp/E7AITD/rWLKAh4hUAlpanpRDCIGxgw7icqdQwIVFL27XN3P1EpYUCgkjJ9urv7iUoJAwKVlOZmoKYm+b6aGuN+olLHgEAlJRwG1q0DQiFAxPi5bh0HlIkAZhlRCQqHGQCIzPAKgYiIADAgEBFRDAMCEREBYEAgIqIYBoSgYsEdIsozZhkFEQvuEJEPeIUQRCy4Q0Q+YEAIIhbcISIfMCAEEQvuEJEPGBCCiAV3iMgHDAhBxII7ROQDZhkFFQvuEFGe8QqBiIgAMCAQEVEMAwIREQFgQKACxwofRLnjW0AQkTNEpFVE9ojIayJyh19tocIUr/DR2QmonqjwwaBAlBk/rxCOA/imqn4WwFwAXxORc31sDxUYVvggyi3fAoKqvqeqr8T+3QdgD4BpfrWHCg8rfBDllqiq322AiNQBeBHADFX9KOWxRgCNADBlypSGlpaWkcf6+/sxduzY/DXUI9yPzESjwODg6PsrK4H6+szek3+L4CiGfQCCsR/z589vV9U5aZ+oqr7eAIwF0A7gy+me29DQoIlaW1u1GHA/MrNxo2pNjaoxgmDcamqM+zPFv0VwFMM+qAZjPwDsUAfHY1+zjESkAsBTACKq+rSfbaFgM8smYoUPotzyrXSFiAiAHwPYo6r/5Fc7KPjSrRfEAECUG35eIfxPAEsAXC4iO2O3RT62hwKK2URE+eFnltFLqiqqOlNVz4/dnvOrPRRcxZhNxAl1FEScqUyBV4jrBdkd8DmhjoKKAYECr9DWC0p3wA9aF1gkGkHdA3Uou6cMdQ/UIRJlZCpVDAgUeIWWTZTugB+kLrBINILGzY3o7O2EQtHZ24nGzY0MCiWKAYEKQjgMdHQAw8PGz6AGAyD9AT9IXWBNW5swcCw5eg0cG0DTVo7YlyIGBKIcS3fAD1IX2L5e8+hldT8VNwYEohxLd8APShdYJAKU9ZtHr+kTAjxiT55hQCCKiUSM+kjZpoI6OeD73QUWH/ge+vdmYDA5etVU1KB5QUBH7MlTDAhUULzK348fIAcHR2cGudlmPGNnyd4yYFUdNuyKBHLMY2TgOxoGNq8DDocAFZT3h7DumnUI1weswZQXvpWuIHIrXQmLbFhlBt1xB3DkiLNtxjN24oO08YwdAIE7wCYNfEfDxg3AsADhe/1pE/mPVwhUMLzM37fKDOrudr7NQsrYCVKmEwUHAwIVDC/z990eCM22WUgZO0HKdKLgYECgguHlWa3VAXLSJOfbtMrMCWLGTlAynShYGBCoYHh5Vhs/QFZWJh8gH3zQ+TabFzSjpqJwMnb8znSi4GFAoILh9VltOGwsvZl4gHSzzXB9GOuuWYfQhBAEgtAE/zJ2enpYTZXcY5YRFRQ/FsRxs81wfdj3jKJIBDh40MiIAnKbjUXFjVcIRFkI4roGTU3GVU4iLihETjAgEDlgduAP6roGQaqmSoWFXUZEaVhNiKuutp6j4GfXDOcYUKZ4hUCUhtWEuO5u8+f7fSbe3GxcySTiHANyggGBKA23B3i/z8TDYSMbinMMyC12GRGlMX36iYydRJMmJdc5AoJzJj5xopE6S+QGrxCI0rCaEPfgg8Ux25drKlMcrxCI0ogf4JuajO6j6dONIBG/v9ACQKJCqtBK3uMVQqaCmIBOnrEq81DoH4NCqtBK3uMVQia8LMxPBaMYPgaFVKGVvMcrhEx4WZif8ir1DL+nx/lri+FjUEgVWsl7DAiZ4FTQomA20/idd4wBYifdP8XwMSi0Cq3kLV8Dgog8LiIHReRVP9vhGqeCFgWzM/w4J2UoiuFjEKQKreQ/v68QngBwlc9tcI/LTRWFdGfy6bp/zD4GIkYwKaQB5nB9GB2rOjC8ehgdqzoYDEqYrwFBVV8E4KLX1gOZpIlwuami4ORM3i5oJH4MAOOjoGr8OyiF7ojcEI1/gv1qgEgdgH9R1RkWjzcCaASAKVOmNLS0tIw81t/fj7Fjx2a+8Z4e45ubWCu4rMz4hk+cmPn7upT1fgREoe2H2Z+/trYfXV0n9qGy0lg0J51oFBgcHH2/09fnmpO/RU8PsH+/0e7KSmDatLx+7NMqtM+TlSDsx/z589tVdU7aJ6qqrzcAdQBedfLchoYGTdTa2qpZCYVUjZO65FsoNPq5Gzca94sYPzduzG7bCbLej4AoxP2I/1kB4097332tIx+Dmhrnf2YR84+SiKfNt5Tub7Fxo7F/iW11s78j77N7o4buD6msEQ3dH9KNu/m9SBWE/QCwQx0cY/0eQ/CX0zSRoBa+p6zFJ5ypAhs2jF5T2UkvYCQyurpoXFAHmHORMhuf5dzZ2wmFjsxyZumLwlXaAcFpmkgxJJxTWmZrKqcTP1cYGhr9WJDzDHKRMstZzsXH77TTnwPYBuAzItIlIl/JawOcZgsVQ8I5ecIqdbW8PNh5BhMviwCr6oDVZcbPeuOs3s0VDWc5Fx+/s4z+QlWnqmqFqtaq6o/z2gCn2ULFkHBOWUlNRlu50vhpVhYbMK4yghoMItEI+uY3Aqd0AqLGz2saUdEQcXVFw1nOxae0u4wA66pliTjvoKSZDSE9/LB1MACCfa7QtLUJg5pyWVM5gPHXN7kKYpzlXHwYEJxwM++g0Mtf0ih2M5rNBP1cwapLp+e4u64eznIuPgwITjm5knCSjcSAUXDcDBUFfY7iyl+vhMJ87lEmXT1ezXKORCOIHoxy0Z48Y0DIpXTZSCtXAkuWjA4YbkpsUt457f4JhZxnJ/lh5a9X4uEdD5s+FqSunng66+DQINNZ84wBIZfsspEiEeCRR07UNogbGDCmi1JgmQ0hpQpiN1HqxeijO9ZZPjdIXT1MZ/UPA0Iu2WUjNTWNDgZxg4PsRgowsyGkr3412KWsenpG914Oq8lkiZigBAOA6ax+YkDIVuJpWH+/MdU1UfzUMV1HNGdBB1rqENLatcafdfp040/b1BSsoaL9+00GwofLTZ9bLub3+4XprP5hQMhG6iByd7fxc9Kk0aeObvIQOQs68OzyB/yqdBKJRlD3QB3K7inD4KnRkclmI3Y0wmw8ubGh0duGucR0Vv9wTeVsmA0iHzsGjB0LHDqUfH9zc/ICvEByveRUnAUdaOnyB6we86pbKT4QO9L3Xj4IXBM70EdjG/3XtTh5HHD0vHUY0iGUSzkaGxqx9otrvWlUhuLdVz17eiAQTJ8wHc0LmgPVrVWsGBCy4aSkRSRiHAn27TNqC1dXGx2806cbQcLqSiDIM5soo2omXsZ4s4FYVA4AC5pGAkJNDfDotWsRDgcrAJgJ14fR1t2G4f89nP7JlDPsMkrlpvPXbhA5EgEmTwZuvjm5S+nIEaOsZjw/sbl5dKnMIKasUBK7P70flU4sB1wn7HM18J3Y7cT8/9LDgJDIbeevVUmLRYuM13V3j37NwACwdOmJ9wyHjW9rkFNWaBS7aiZ+VDqxGnANnTLdcfXWdOWsGSyKHwNCIrdlrq1KWjz3nH2tg6Gh5EAzcWL6WdDkmcSLwmg0+1VUzR5butT4GFldeGablZSLgVi7/H8Gi9LAgJAoXcew2bfWrKSFk85iZhIFQupF4eCgs4ygxKGh+HBQYhxP/Fg0NwPr11tfeOYiKym1rlBleaXryWZ2+f/ZBAsqHAwIidKNCTj91jrtLGYmke8yWfvI7QE83TZytf5SYl2h+k/Uu87Kscv/zzRYUGFhQEhk1/lr9a2NjwekTlCrqEi/PWYS+S6TbCG3B/B02wjK+kt23U6ZBgsqLAwIiew6hq2+nUNDwLJlwG23JWcTiZyYoDZpkvUMZvJVJhlBbg/g6baR76wkq/5+u3LWmQYLKiwMCKmsylzbfTuPHTM6nxMNDhoT1IaHjUlqjz/OTKIAyiQjaOJE8/utPiLptpGPrKR4EJB7BEueXmLZ329VzjrTYEGFhQHBKSclL1MlnjI6WU+B8i71orCy0j5WRyLARx+Nvr+y0voAnm59JTfrL2UicdAXwKj1EJz292cSLKiwcKayU/Fv59KlRjeRExwjKAjxVFEAaGsD5s2zfm5Tk3FBmGrcOPsDeOI2rB6Pv3+8WF7i/dkwncWcItv+/nB9mAGgCPAKwY1w2MgfTL1SqKgYPUYgYowpsJx1QbKaF2A1ThBf4yjT+QReFcSLRIDOw+kP9uzvJyBNQBCRc0RkgYiMTbn/Km+bFWBm1/c/+cmJMQIguWgdy1kXHLOD85IlxoJ3ucpMTmWVuXTzzZmfU8Tbg177g32++vtTB7NX/nolJ7MFjGVAEJFvAHgWwO0AXhWR6xIe/oHXDQs0s/GA+H2hkPmqaJyEVjDMDs6qxoJ3ixa5z0y2+tMnXk10dlq3J9NzipH2bG0GBpMbLRAAyFt/v9nktYd3PMzJbAFjd4XwlwAaVPVLAOYB+I6I3BF7TLxuWMEKSlI5ZczqT6VqVCVxm5lsdn/q1UQ6mZxTjGw3GgY2rwMOhwAV4HAIG768AbpakwaHveRkHIOT2fxnN6hcrqr9AKCqHSIyD8AmEQmBAcHa9Onmp3scYC4YVn9CwDjIWg0Qu/nTm11NpOP2nCKpPdHwSBnsUAgI3+/uvbLldNCak9n8ZXeF8L6InB//JRYcrgYwGUC91w0rWH6UuqScam42zv7N2MV1sz99ZaUxcd3p4LQdt+cUTj+K+Vju0+mgNQe3/WUXEG4B8H7iHap6XFVvAXCpp60qZF4nlZPnwmFgxYrRQSFdXE/900+alLyyauJYgNXB3WxSO2Aksrk9p3DyUXQ6EJ5tNVOzyWupOJnNf5YBQVW7VPV9i8f+MxcbF5GrRORNEdkrIt/KxXsGAiehFby1a411jNzG9cQ//dixo+csxMcCrM7egdGT3gFg/Hj7bacesHuO9Ixqj9lH0clAeC6qmZpNXvvqnK9yMlvA+DYxTUTKATwE4AoAXQC2i8ivVPV1v9qUM+lqI1NBSDeZLB27QebEiWidnUB5uXEgthpXiM9zSBWJRnDHv96B7iMnFmPq7O1EZ28nItFI2gOsk4Fwu2qmbg7gnLzmTiQaQdPWJuzr3Ze3daX9nJj2OQB7VfUPqjoIoAXAdWle4w83naxezTCigpOuaF18BdWamvST383eKxKN4LZnGpOCQdywDjvK2HFSWI/VTPPPrzUmRC1y3kTkLABTUruHROR/ATigqr/PasMiNwK4SlWXx35fAuAiVf16yvMaATQCwJQpUxpaWlpGHuvv78fYsUlz5nKvp8c4qA8nLPZdVmb0IZhVOYtGza/5KyuBevOx+LzsRx4Uw37kch+cfHSsPi6JrD5uu96P4riav7j2pFp0fdyFyvJKTBs3DROrzSvyOWrjwSgGh0Zvp7K8EvWf8C6/pBg+T0Bm+5Hr//P58+e3q+qctE9UVdMbgH8BMNPk/jkANlu9zukNwGIAjyX8vgTAD+1e09DQoIlaW1vVc6GQqnGun3wLhcyfL2L+fBHLTdjux8aNxrZEjJ8bN2a+Lx7Ly9/DY7neh3R/PquPS/wjY/cnx2pRrIHp7b6f3Tfy75rmGt242/pzk66NG3dv1JrmmqT3T/eeuVAMnyfVzPZD1pj/bWWN9XHEDoAd6uC4bNdlVKequ00CyA4Ada7Ck7kuAGck/F4L4EAO3je37DpZzbqS3NZGtsPup4KXblDX6mMRCjnISUhTkiIu3YSvdG1kNdP882uNCbuAUGXzWHUOtr0dwKdF5EwRqQRwE4Bf5eB9c8vqGztx4uiD9bJlwOHDo59rVxvZTq7WVqTAymbayqSdo0tSwLwHOCfVTM1KX5M3/Fpjwi4gbBeRv0y9U0S+AqA92w2r6nEAXwfw7wD2APiFqr6W7fvmnF1+YOrB+tgx89HBdLWRrbAMRtHLZtrKg8vDqPj35JIUODLJ9Lmc8FVY/Loqs0s7XQXgGREJ40QAmAOgEsD1udi4qj4H4LlcvJdnUgvVx9NIlyxx/h5WOYPpsAxGScg0vdV4TRhNTeGRj+aiuyJY/8fGpOdxwldh8iNN125i2geqegmAewB0xG73qOrFajFhrWiZdbK6OSgnPtdNCivLYFAaqR/NtV81ziwryyvZ30+u2ZW/rhKRVQBuADAI4GFVfSFvLQs6N0tqxhfKWbnSfJDY6gqCZTDIgt15Rbg+jPpP1Oe8vz9xNvTkf5iMyf8wuaDWMsi2/EYpsOsyWg/gGIDfAPgzAJ+F0Y1EwOiupLIy+9lFnZ1GQX2ztRL277ffDgMAJYgnn8WHsOLnFYB3H5X4RKn4jOXUmdGNm40GuAk++ZyJm9r+TNtc7OwGlc9V1ZtV9VEAN4IF7UZLvF5PnNljxarwfbqZSUQJ/Eg+S7eegdu1DPI9E9eu/AadYBcQRspyxTKCyE42A71m5S2JLOQj+Sy1S6rTQdqqm9TWfB+gWX7DGbuAMEtEPord+gDMjP9bRD7KVwMLhpsxhUQ1NcC0ablvDxUtJ/WHsmE2H1IcTIJzk9qa7wO0XxO9Co1dllG5qo6P3cap6piEf4/PZyMLglkxfDuJg8Tx2c35WKmECobVx8Hr5DPTNaW3NEOOW5/wuE1tzfcB2q+JXoXGz2qnxSdxTOHQIeugYFaXwOy0bNkyYPJkBogSZFe1xOvks87xEWBVHbC6zPhZHwGiYeizJyZKTaqehEnVkzJObV306UWm95818ays2m6VSeT1RK9iyWDybT2EvPJ6fQKr93/wweR0EMD6VM7stOzYMWO5LSA/qSQUGHYDx/HEMy8+BpFoBHJdI3RMbOOndALXGJ+70EdhdKzKzUafe9t8PuoL77zgaB0HM+kyibya6FVMGUzFf4XgdYE4s/dfsgRYuPDEt7q83Hiu1alcvAZxOqxjVDL8qlrStLXpRDCIqxyALGzK6XxIq7EChWY8sOxXJlExZTAVf0DwOkfPtMNVga1bTxzkh4ZOXBmkBoNIxFkwiGMdo5KQbuA43XBTT09mw1GWB+oJ+3J6RWI3VpDpwLJfmUTFlMFU/AHB61Mtp+8zMADcccfo+5uanM1hiGMdo5JgN3Cc7qI3fo6RyUWx1YE6lOPB3uYFzRCIqzak41cmUTFlMBV/QPA6R8/N+3R3j/5W2gWU1PkJrGNUMuwGjtNd9JqdYzi9KM5XNk64PozLz7x81P1ut5U4mNs/2I/K8uTvTD4yiYopg6n4A0ImOXpuC9CJ+ZmOqdRvpd0KKY8/zjpGNEq6i95sLorzVXY5Eo1gW9e2pPsEgqWzljreVups5+4j3VDVrLKfMlFMCwgVf5aRVflqqwOr20Ix4TDwn/9pXqfITOq3srkZOHgw+b7E8QYGgJJk9zFMVxXdzUWxeYKc92WXzQZiFTqSfRSJRtBzsAeX33O5ZZ0js/c4NnwMYyvH4tD/OeRp+1P5UaraC8V/hQCkXyMwUbrrcbOrh7VrgQ0bks/mrRbVTv1WhsPG83klQAnsPobpLnqbm42Pp9XjcX6u0Go3EBs/8x8cGrStc1RMg7lBURoBwY10ayjbzRZKDDqPPGI/KpgYVADnAYtKgt3HMN3ENKfnGLlKwMtkUtbEavO1x6dPmO44jbOYBnODggEhld31tptvkNW3FhgdVDo7OQuZkqTr9kl30TtxYvpzDNtzH4cH+UyqlkaiEfQN9pk+1j/Yj85e8zTs1DP/YhrMDQoGhFR21+NuR+vMvrVmQWV4mBPOKEk+FsuzCjoTL3N+kM9kUlbT1iYMDpmXfO8+0u04HbWYBnODggEhld31eC5SWP2agkoFJR+L5VkFHSx0fpDPpB8/XR+/QkcFBasz/3B9GB2rOnK+OlypYkAwY3U9novTNq/nRRA5ZBV0eo47P8hn0o/vpI9foVwX2gcMCG7k4rTNLKiUlXHCGSXJVwaQ2bmPm4N8Jv34Zq9JFZoQ8mRd6GwUS0VTOwwIbrlJYbV6fWpQCYWYWVSkMl3iwo9lMuPcHOQz6cdPfA0Ax91Dfsr3kp9+YUDwQ2pQmWiegkeFLZuz/FwONbkNSm4P8tn04wsEE6snms4u7jnSE5gzcqvB86XPLC2qoFD8M5WJfJJuTQM76WYjO+V24n2clzNvU9cP6D7SjZqKGmz48oaRbUaiERzsPTiSgur3GgNWA+FDOlSwax+Y4RWCF7gUJiG7s/xcpZ362fVkxUmqatPWJgzrsO1z8sluILxQ1z4ww4CQa37WA6BAySahLFdpp0HMcnaSqmo1Oc3qfq+lGwgvlnIZDAi5FsRTMvJFtmf52eYvAMHMcnaSxVQu5abPsbrfa/FxFavtF0u5DF8CgogsFpHXRGRYROb40QbPWK1+xolnJScfk8vSyTQo5SLF0uo9rLKYFn160cjzh3TI9D2t7s+HcH0Y669fX9TlMvwaVH4VwJcBPOrT9r0RiRjffLMy2GVlJ4rgUcnwu4K52+rvQG4WjXfyHk1bm7Cvdx+mT5iORZ9ehPW71o8aW0gVT1X1i1nbzUpzFypfAoKq7gEAcbOwTCFoarJeE2FoyFl6B1GOuQ1KdoO+qQe+SDRienBM9x6pWUx1D9SlDQZBORMvlrUPzIg6WdTFq42LtAG4U1V32DynEUAjAEyZMqWhpaVl5LH+/n6MtVp3wA/t7emfU1kJ1Ncn3RW4/chQMexHMewDkN1+tL9n/TlumNow8u+eIz3o7O1MygYqkzKEJoTwzuF3HL2H3TZrT6pF18ddAIAxZWNwxvgzLMtmB1kQPlPz589vV9W03fOeBQQR2QLgdJOHmlT12dhz2pAmICSaM2eO7thx4qltbW2YN29e9o11I77EVGcnUF5unPmHQsZ1ePx+OyKjFrz1ZT88UAz7UQz7AGS3H3UP1Jlm84QmhNCxqsPR8wDzjKDU97B7r/vOvg93vnWn7esKQRA+UyLiKCB4NqisqgtVdYbJ7Vmvtum5xJRSwAgGwInU0kWLRo/gpeKsZAo4p6Ur7NJH3dY4KpW0zqBj2qkbZimlcQMDwHPPnUgrsdLXxzkJFGhOS1fYpY9mUv6iFNI6g86XQWURuR7ADwGcBuDXIrJTVa/0oy2upEsdja9vGB/BmzwZ6O5Ofs7goLPaBUQ+cjJw2rygOSmTCEi+CnA7+Bp/rt17krd8uUJQ1WdUtVZVT1LVKQURDID0s3lSH+/pMX8e5yRQEfBixbLUSqhcCyG/WNzOjebm5Ephicxm++SqQhlRQHmRghl/z7a2NnT8RUdO35vscQzBjcSpp4CRZQRYT0HNx8K4REQ5wisEt9zM8slkmigRkU8YELzmd+0CIiKH2GVEREQAGBCIiCiGASFouNoaEfmEYwhBkukCuEREOcArBC9kepbP1daIyEcMCLlmt6ZyukARxAVwiahksMso16zO8u+4AzhyxLw7aNo04ydnNhORj3iFkGtWZ/Pd3em7gzizmYh8xICQa27P5hMDSBBWZSeiksUuo1wzK4BXUwNUV48uhQ2MDiCc2Zzk2LFj6OrqwtGjR/OyvQkTJmDPnj152ZaXrPajqqoKtbW1qKio8KFVFHQMCLlmVb8IMA8U7A6y1dXVhXHjxqGurg4i4vn2+vr6MG7cOM+34zWz/VBVdHd3o6urC2eeeaZPLaMgY0Dwgt1Zvlmhu7a2vDavkBw9ejRvwaDYiQgmTZqEDz/80O+mUEAxIOQTu4MywmCQO/y/JDscVCYiIgAMCEQ59cQTT+DAgQN+N4MoIwwIVFx8Lg7IgECFjAEB8P0gQjliVzYkC3/605/wxS9+EbNmzcKMGTPw5JNPor29HZdddhkaGhpw5ZVX4r333sOmTZuwY8cOhMNhnH/++Thy5Ai2bt2K2bNno76+Hrfddhs+/vhjAMC3vvUtnHvuuZg5cybuvPNOAMDmzZtx0UUXYfbs2Vi4cCE++OCDrP9LiFxR1YK5NTQ0aKLW1lbN2saNqjU1qsYhxLjV1Bj350lO9iMAvNiP119/3fmTQ6Hkv2P8Fgo5fouPPvpo1H2bNm3S5cuXj/x++PBhvfjii/XgwYOqqtrS0qLLli1TVdXLLrtMt2/frqqqR44c0draWn3zzTdVVXXJkiV6//33a3d3t5599tk6PDysqqp//OMfVVW1p6dn5L4f/ehH+td//dfO993BfsS5+j/1Eb8XuQNghzo4xvIKwYsKo7zi8IdHxQHr6+uxZcsW3HXXXfjNb36Dd999F6+++iquuOIKnH/++fjbv/1bdHV1jXrdm2++iTPPPBNnn302AGDp0qV48cUXMX78eFRVVWH58uV4+umnURMrV9LV1YUrr7wS9fX1uPfee/Haa69l1W4itxgQcn0Q8ajbghywKhuSZXHAs88+G+3t7aivr8fdd9+Np556CufIk6eWAAAPx0lEQVSddx527tyJnTt3IhqN4vnnnx/1OuPEbLQxY8bg5Zdfxg033IBf/vKXuOqqqwAAt99+O77+9a8jGo3i0UcfzdvsbKI4BoRcH0S4poF/PCoOeODAAdTU1ODmm2/GnXfeid/97nf48MMPsW3bNgBGeY342fy4cePQ19cHADjnnHPQ0dGBvXv3AgA2bNiAyy67DP39/ejt7cWiRYvwwAMPYOfOnQCA3t5eTItVvl2/fn1WbS41kWgEdQ/UoeyeMtQ9UIdIlCdgmeDENKvaQ5keRLimgX+syoZkORkwGo3ib/7mb1BWVoaKigo8/PDDGDNmDL7xjW+gt7cXx48fx6pVq3Deeefh1ltvxYoVK1BdXY1t27bhJz/5CRYvXozjx4/jwgsvxIoVK9DT04PrrrsOR48ehari/vvvBwCsWbMGixcvxrRp0zB37ly888472f6PlIRINILGzY0YOGZ8hzt7O9G42SgtH67nRFBXnAw0BOXmyaCyqjGAHAqpihg/sxlQzmBgMwiDTrng+6ByDtgNxhaSUhpUDt0fUqzBqFvo/pCn7XMqCN9vBHlQWUTuFZE3RGS3iDwjIqf40Y4R4TDQ0QEMDxs/szmj5JoGRHm1r9f86tvqfrLm1xjCfwCYoaozAbwF4G6f2pF7XNOAKK+mTzAf77O6n6z5EhBU9XlVPR779bcAav1oh2dyecVBRLaaFzSjpiL5qrymogbNC3hV7paoRWpc3hogshnAk6q60eLxRgCNADBlypSGlpaWkcf6+/sxduzYvLTTS9wPaxMmTMBZZ52V0/e0MzQ0hPLy8rxtzyt2+7F371709vbmuUXuufk89Rzpwf6+/RgcGkRleSWmjZuGidUTPW6hM0H4fs+fP79dVeeke55nWUYisgXA6SYPNanqs7HnNAE4DsAyR0xV1wFYBwBz5szRefPmjTzW1taGxN8LFffD2p49e/K6YE0xL5ATV1VVhdmzZ+e5Re7xe5F/ngUEVV1o97iILAVwNYAF6vdlChER+ZZldBWAuwBcq6oD6Z7vO5aioBz67ne/iy1btrh+XVtbG66++moPWkRk8CvL6P8CGAfgP0Rkp4g84lM70mMpioISlNitqhgeHjZ97Hvf+x4WLrS9gM6J48ePp38SUQK/sozOUtUzVPX82G2FH+1whKUoCoYXsfuuu+7C2rVrR35fs2YN/vEf/xH33nsvLrzwQsycOROrV68GAHR0dOCzn/0sVq5ciQsuuADvvvsubr31VsyYMQP19fUjM5JvvfVWbNq0CQCwfft2XHLJJZg1axY+97nPoa+vD0ePHsWyZctQX1+P2bNno7W1dVS7enp68KUvfQkzZ87E3LlzsXv37pH2NTY24rrrrsMtt9yS+Y5TSWIto3RYiqJgeBG7b7rpJjz55JMjv//iF7/Aaaedhrfffhsvv/wydu7cifb2drz44osAjAqnt9xyC/77v/8bhw4dwv79+/Hqq68iGo1i2bJlSe89ODiIP//zP8eDDz6IXbt2YcuWLaiursZDDz0EwCiZ8fOf/xxLly4dVehu9erVmD17Nnbv3o0f/OAHSQf/9vZ2tLS04Gc/+1nmO04liQEhHY8qaFLueRG7Z8+ejYMHD+LAgQPYtWsXTj31VOzevRvPP/88Zs+ejQsuuABvvPEG3n77bQBAKBTC3LlzAQCf+tSn8Ic//AG33347/u3f/g3jx49Peu8333wTU6dOxYUXXggAGD9+PMaMGYOXXnoJS5YsAWAUyAuFQnjrrbeSXpv4nMsvvxzd3d0jqaTXXnstqqurM99pKlkMCOmwFEXB8Cp233jjjdi0aROefPJJ3HTTTVBV3H333SPlr/fu3YuvfOUrAICTTz555HWnnnoqdu3ahXnz5uGhhx7C8uXLk95XVSEio7bnJOnO7Dnx90psA5EbDAjpsBRFwfAqdt90001oaWnBpk2bcOONN+LKK6/E448/jv7+fgDA/v37cfDgwVGvO3ToEIaHh3HDDTfg+9//Pl555ZWkx8855xwcOHAA27dvB2DMHTh+/DguvfRSRGIDH2+99Rb27duHz3zmM0mvTXxOW1sbJk+ePOoKhMgtlr92IhxmACgAHlW/xnnnnYe+vj5MmzYNU6dOxdSpU7Fnzx5cfPHFAICxY8di48aNo2YG79+/H8uWLRvJNvq7v/u7pMcrKyvx5JNP4vbbb8eRI0dQXV2NLVu2YOXKlVixYgXq6+sxZswYPPHEEzjppJOSXrtmzRosW7YMM2fORE1NDddPoNxwUhI1KDfPyl/7jPthjeWvM1NK5a+DLgj7gSCXvyYiouBhQCAiIgAMCEREFMOAQEREABgQiIgohgGBiEpCJBpB3QN1KLunDHUP1CESZYHKVAwIRC4dOHAAN954o+vXLV++HK+//rrtcx555BH89Kc/zbRpZCESjaBxcyM6ezuhUHT2dqJxcyODQgoGBCoq+TgL/OQnPzlSrTRRunLTjz32GM4991zb56xYsYJVSj3QtLUJA8eSKx8OHBtA01ZWLU7EgEBFw4uzQKvy1zNmzAAAPPHEE1i8eDGuueYafOELX8Dw8DBWrlyJ8847D1dffTUWLVo0EjzmzZuHHTt2ADBmNzc1NWHWrFmYO3cuPvjgg5H3v++++wAYax8vXLgQs2bNwgUXXIDf//736O/vx4IFC3DBBRegvr4ezz77bMb7VqzMTgr29ZpXOLS6v1QxIFDR8OIs0Kz8dbw6ady2bduwfv16vPDCC3j66afR0dGBaDSKxx57DNu2bTN93z/96U+YO3cudu3ahUsvvRQ/+tGPRj0nHA7ja1/7Gnbt2oX/+q//wtSpU1FVVYVnnnkGr7zyClpbW/HNb37TUTG8UmF1UjCxeqLp86dPYNXiRAwIuRKUpbpKmBdngWblr6enlE+94oorMHGiccB56aWXsHjxYpSVleH000/H/PnzTd+3srJyZDnMhoYGdHR0JD3e19eH/fv34/rrrwcAVFVVoaamBqqKb3/725g5cyYWLlyI/fv3j1xdkPVJAQDUVCRXPqypqEHzAlYtTsSAkAtcZjMQrM72sj0LTC1/nSqx3LTTs/WKioqRctXl5eWjxh+s3icSieDDDz9Ee3s7du7ciSlTpoxaPKeUWQX/niM9WHfNOoQmhCAQhCaEsO6adQjXs2hlIgaEXOAym4HQvKDZk7PA1PLXdj7/+c/jqaeewvDwMD744AO0tbVltM3x48ejtrYWv/zlLwEAH3/8MQYGBtDb24tPfOITqKioQGtrKzo7OzN6/2Jld1IQrg+jY1UHhlcPo2NVB4OBCQaEXOAym4EQrg97chaYWv7azg033IDa2lrMmDEDf/VXf4WLLroIEyZMyGi7GzZswD//8z9j5syZuOSSS/D+++8jHA5jx44dmDNnDiKRCM4555yM3rtYeXVSUDKclEQNyi2w5a9DIVWjsyj5Fgo5enlg9iNLLH9t6OvrU1XVQ4cO6ac+9Sl97733sn5Pt0q5/PXG3Rs1dH9IZY1o6P6Qbty9MbcNcykI3284LH/NBXJyobnZGDNI7DbiMpsl6+qrr8bhw4cxODiI73znOzj99NP9blJJCdeH2R2UIQaEXPBqqS4qSJmOGxD5jQEhV7jMpmfUYjF6ck85Z4FscFCZAq2qqgrd3d08kOWAqqK7uxtVVVV+N4UCilcIFGi1tbXo6urChx9+mJftHT16tCgOmFb7UVVVhdraWh9aRIWAAYECraKiAmeeeWbettfW1obZs2fnbXteKZb9oPxilxEREQFgQCAiohgGBCIiAgBIIWVviMiHABKLt0wGcMin5uQS9yM4imEfgOLYj2LYByAY+xFS1dPSPamgAkIqEdmhqnP8bke2uB/BUQz7ABTHfhTDPgCFtR/sMiIiIgAMCEREFFPoAWGd3w3IEe5HcBTDPgDFsR/FsA9AAe1HQY8hEBFR7hT6FQIREeUIAwIREQEogoAgIt8Xkd0islNEnheRT/rdpkyIyL0i8kZsX54RkVP8bpNbIrJYRF4TkWERKYg0u0QicpWIvCkie0XkW363JxMi8riIHBSRV/1uS6ZE5AwRaRWRPbHP0x1+tykTIlIlIi+LyK7Yftzjd5vSKfgxBBEZr6ofxf79DQDnquoKn5vlmoh8AcALqnpcRP4eAFT1Lp+b5YqIfBbAMIBHAdypqjt8bpJjIlIO4C0AVwDoArAdwF+o6uu+NswlEbkUQD+An6rqDL/bkwkRmQpgqqq+IiLjALQD+FIB/i0EwMmq2i8iFQBeAnCHqv7W56ZZKvgrhHgwiDkZQEFGOFV9XlWPx379LYCCq1GsqntU9U2/25GhzwHYq6p/UNVBAC0ArvO5Ta6p6osAevxuRzZU9T1VfSX27z4AewBM87dV7sWWM+6P/VoRuwX6+FTwAQEARKRZRN4FEAbwXb/bkwO3AfhXvxtRYqYBeDfh9y4U4EGo2IhIHYDZAH7nb0syIyLlIrITwEEA/6Gqgd6PgggIIrJFRF41uV0HAKrapKpnAIgA+Lq/rbWWbj9iz2kCcBzGvgSOk30oUGZrdAb6bK7YichYAE8BWJXSE1AwVHVIVc+HccX/OREJdDdeQSyQo6oLHT71ZwB+DWC1h83JWLr9EJGlAK4GsEADOrjj4m9RaLoAnJHwey2AAz61peTF+tyfAhBR1af9bk+2VPWwiLQBuApAYAf8C+IKwY6IfDrh12sBvOFXW7IhIlcBuAvAtao64Hd7StB2AJ8WkTNFpBLATQB+5XObSlJsMPbHAPao6j/53Z5Michp8WxBEakGsBABPz4VQ5bRUwA+AyO7pRPAClXd72+r3BORvQBOAtAdu+u3hZYtJSLXA/ghgNMAHAawU1Wv9LdVzonIIgAPACgH8LiqNvvcJNdE5OcA5sEoufwBgNWq+mNfG+WSiHwewG8ARGF8rwHg26r6nH+tck9EZgJYD+PzVAbgF6r6PX9bZa/gAwIREeVGwXcZERFRbjAgEBERAAYEIiKKYUAgIiIADAhERBTDgEBkQ0SGYpV0XxWR/yciNbH7TxeRFhH5vYi8LiLPicjZJq8v+OqjVDoYEIjsHVHV82OVQwcBrIhNnHoGQJuq/g9VPRfAtwFMMXn9EzBmpxIFXkGUriAKiN8AmAlgPoBjqvpI/AFV3Wn2AlV9MVagjSjweIVA5ICIjAHwZzBmz86AUaOfqKgwIBDZq46VL94BYB+MGjtERYldRkT2jsTKF48QkdcA3OhTe4g8wysEIvdeAHCSiPxl/A4RuVBELvOxTURZY0Agcim2VsX1AK6IpZ2+BmANTNZPiFUf3QbgMyLSJSJfyWtjiVxgtVMiIgLAKwQiIophQCAiIgAMCEREFMOAQEREABgQiIgohgGBiIgAMCAQEVHM/wePf4sNfXU7WwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for label, colour in zip(classes, colours):\n",
    "    print(label, colour)\n",
    "    plt.scatter(projections[label_names == label, 0],\n",
    "                projections[label_names == label, 1],\n",
    "                label = label,\n",
    "                c = colour)\n",
    "plt.xlabel('PC 1')\n",
    "plt.ylabel('PC 2')\n",
    "plt.legend(loc='lower center')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Related materials\n",
    "* [Wikipedia](https://en.wikipedia.org/wiki/Principal_component_analysis)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nav_menu": {},
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}