Principal Component Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![tinhatben](tinhatben_svg.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis\n",
    "\n",
    "This Jupyter notebook has been written to partner with the [Principal components blog post](https://tinhatben.com/2016/04/13/principal-component-analysis/)\n",
    "\n",
    "Copyright (c) 2016, tinhatben\n",
    "All rights reserved.\n",
    "\n",
    "Redistribution and use in source and binary forms, with or without\n",
    "modification, are permitted provided that the following conditions are met:\n",
    "\n",
    "* Redistributions of source code must retain the above copyright notice, this\n",
    "  list of conditions and the following disclaimer.\n",
    "\n",
    "* Redistributions in binary form must reproduce the above copyright notice,\n",
    "  this list of conditions and the following disclaimer in the documentation\n",
    "  and/or other materials provided with the distribution.\n",
    "\n",
    "* Neither the name of tinhatben nor the names of its\n",
    "  contributors may be used to endorse or promote products derived from\n",
    "  this software without specific prior written permission.\n",
    "\n",
    "THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\"\n",
    "AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\n",
    "IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE\n",
    "DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE\n",
    "FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL\n",
    "DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR\n",
    "SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER\n",
    "CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,\n",
    "OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE\n",
    "OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------------------------------------------------------------------------------------------------------------\n",
    "-----------------------------------------------------------------------------------------------------------------------\n",
    "# Components of PCA\n",
    "\n",
    "## Mean\n",
    "\n",
    "In this example we will show an implementation of an equation to calculate the mean, as well as using Python and numpy's interal functions.  We recall from the blog post:\n",
    "\n",
    "$$ \\large{\\mu = \\frac{1}{n}\\sum_{k=1}^n x_k }$$\n",
    "\n",
    "The numpy documentation for numpy.mean can be found [here](https://docs.scipy.org/doc/numpy/reference/generated/numpy.mean.html)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Very basic data: [1 2 3 4 5]\n",
      "Very basic data 2D:\n",
      " [[ 1  2  3  4  5]\n",
      " [ 6  7  8  9 10]]\n",
      "\n",
      "----Calculating By Hand----\n",
      "\n",
      "The mean of very_basic data is : 3.00\n",
      "The mean of very_basic data - 2D is : 3.00 and 8.00\n",
      "\n",
      "----Calculating Using NumPy----\n",
      "\n",
      "The mean of very_basic data is (calculated using NumPy): 3.00\n"
     ]
    }
   ],
   "source": [
    "# Calculating the mean by hand on some very simple data\n",
    "very_basic_data = np.array([1,2,3,4,5])\n",
    "very_basic_data_2d = np.array([[1,2,3,4,5],[6,7,8,9,10]])\n",
    "\n",
    "print(\"Very basic data: %s\" % very_basic_data.__str__())\n",
    "print(\"Very basic data 2D:\\n %s\\n\" % very_basic_data_2d.__str__())\n",
    "\n",
    "n = very_basic_data.shape[0]\n",
    "n_2d = very_basic_data_2d.shape[1]\n",
    "\n",
    "mu = (np.sum(very_basic_data)) / n\n",
    "\n",
    "mu_2d = (np.sum(very_basic_data_2d, axis=1)) / n_2d\n",
    "\n",
    "print(\"----Calculating By Hand----\")\n",
    "\n",
    "print(\"\\nThe mean of very_basic data is : %0.2f\" % mu )\n",
    "print(\"The mean of very_basic data - 2D is : %0.2f and %0.2f\\n\" % (mu_2d[0], mu_2d[1]) )\n",
    "\n",
    "print(\"----Calculating Using NumPy----\")\n",
    "\n",
    "# Calculating using numpy\n",
    "\n",
    "mu = np.mean(very_basic_data)\n",
    "mu_2d = np.mean(very_basic_data_2d, axis=1)\n",
    "print(\"\\nThe mean of very_basic data is (calculated using NumPy): %0.2f\" % mu )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standard Deviation\n",
    "Recall that he standard deviation is:\n",
    "\n",
    "$$ \\large{\\sigma = \\sqrt{\\frac{\\sum_{k=1}^n (x_k - \\mu)^2}{n - 1}}}$$\n",
    "\n",
    "The numpy documentation to the standard deviation  (numpy.std) can be found [here](https://docs.scipy.org/doc/numpy/reference/generated/numpy.std.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "----------------------------------------1D Data----------------------------------------\n",
      "The very basic data set with some random noise: [ 5 10 15 20 25]\n",
      "The new mean (mu) = 15.00\n",
      "Variance (sigma^2) = 62.50\n",
      "Standard deviation (sigma) = 7.91\n",
      "Standard deviation using NumPy (sigma)  = 7.91\n",
      "\n",
      "----------------------------------------2D Data----------------------------------------\n",
      "The very basic data set (2D) with some random noise: [[ 4  8 12 16 20]\n",
      " [ 6  7  8  9 10]]\n",
      "The new 2D mean (mu) = 12.00 and 8.00\n",
      "Variance (sigma^2) = 40.00 and 2.50\n",
      "Standard deviation (sigma) = 6.32 and 1.58\n",
      "Standard deviation using NumPy (sigma)  = 6.32 and 1.58\n",
      "\n",
      "\n",
      "**Note: The variance calculated for the 2D data sets, considers each row or array as an independant data set\n"
     ]
    }
   ],
   "source": [
    "# Add some noise to very_basic_data to ensure \n",
    "from random import randint\n",
    "\n",
    "very_basic_data *= randint(1,5)\n",
    "very_basic_data_2d[0,:] *= randint(1,5)\n",
    "very_basic_data_2d[1,:] *= randint(1,5)\n",
    "\n",
    "mu = np.mean(very_basic_data)\n",
    "mu_2d = np.mean(very_basic_data_2d, axis=1)\n",
    "\n",
    "print(\"\\n\" + 10 *\"----\" + \"1D Data\" + 10*\"----\")\n",
    "\n",
    "print(\"The very basic data set with some random noise: %s\" % str(very_basic_data))\n",
    "print(\"The new mean (mu) = %0.2f\" % mu)\n",
    "\n",
    "# Using the very basic data\n",
    "sigma_squared = np.sum((very_basic_data - mu) ** 2) / (n - 1)\n",
    "\n",
    "print(\"Variance (sigma^2) = %0.2f\" % sigma_squared)\n",
    "\n",
    "# Standard deviation\n",
    "\n",
    "sigma = np.sqrt(sigma_squared)\n",
    "\n",
    "print(\"Standard deviation (sigma) = %0.2f\" % sigma)\n",
    "\n",
    "# You can also calculate the standard deviation using numpy\n",
    "sigma = np.std(very_basic_data, ddof = 1) # ddof = 1, is used to ensure the denominator used is n - 1, not n\n",
    "print(\"Standard deviation using NumPy (sigma)  = %0.2f\" % sigma)\n",
    "\n",
    "\n",
    "print(\"\\n\" + 10 *\"----\" + \"2D Data\" + 10*\"----\")\n",
    "\n",
    "print(\"The very basic data set (2D) with some random noise: %s\" % str(very_basic_data_2d))\n",
    "print(\"The new 2D mean (mu) = %0.2f and %0.2f\" % (mu_2d[0], mu_2d[1]))\n",
    "\n",
    "# Using the very basic data\n",
    "\n",
    "sigma_squared_2d = np.zeros(2)\n",
    "sigma_squared_2d[0] = np.sum((very_basic_data_2d[0,:] - mu_2d[0]) ** 2) / (n_2d - 1)\n",
    "sigma_squared_2d[1] = np.sum((very_basic_data_2d[1,:] - mu_2d[1]) ** 2) / (n_2d - 1)\n",
    "\n",
    "print(\"Variance (sigma^2) = %0.2f and %0.2f\" % (sigma_squared_2d[0], sigma_squared_2d[1]))\n",
    "\n",
    "# Standard deviation\n",
    "\n",
    "sigma_2d = np.sqrt(sigma_squared_2d)\n",
    "\n",
    "print(\"Standard deviation (sigma) = %0.2f and %0.2f\" % (sigma_2d[0], sigma_2d[1]))\n",
    "\n",
    "# You can also calculate the standard deviation using numpy\n",
    "sigma_2d = np.std(very_basic_data_2d, ddof = 1, axis=1)\n",
    "\n",
    "print(\"Standard deviation using NumPy (sigma)  = %0.2f and %0.2f\" % (sigma_2d[0], sigma_2d[1]))\n",
    "\n",
    "print(\"\\n\\n**Note: The variance calculated for the 2D data sets, considers each row or array as an independant data set\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Covariance\n",
    "\n",
    "The covariance of X and Y is defined as:\n",
    "\n",
    "$$\\large{cov(X,Y) = \\frac{\\sum_{k=1}^n(X_k - \\mu_X)(Y_k - \\mu_Y)}{n-1}}$$\n",
    "\n",
    "At this point we will introduce a new fake data set.  This is a set of sample data I <i>made up</i>, the data consists of two columns and 10 rows.  The first column represents the number of coffees consumed and the second column represents pages of a PhD thesis written."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fc6358b8a90>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAH4CAYAAACxCm+2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYZXV97/v3p2gKlamJCig74uFoNA6IRo0ES8s46wkY\nj/EY4mzUOLREE2+QeJ5uczweNU9w6MhNIupFbzRRr/M8YGlpHJEGBDUmBqQktBEp0IAWsr/3j72q\nLdqu6l3dtWvvVfv9ep799Fprr+G7qqE/9fut31orVYUkSWqniWEXIEmS9p1BLklSixnkkiS1mEEu\nSVKLGeSSJLWYQS5JUosZ5NIISvKcJFcmuTbJEUlOSvLPzfzJw65vNZI8JcnsAPb74yS3X+v9Sm1j\nkEsDkuTUJF9tAuf7ST6c5KQ+ttsE/BXwkKo6rKquBl4GvL6Z/8Cga1+NJN9ofsG4NsnPk1zfnPO1\nSU5vVlvzB1ZU1aFVdela71dqG4NcGoAkLwLOBF4OHAncDngD8Dt9bH40cBDwzSXLjgUuWeMy10RV\n3a35BeMwYBZ4bhOyh1XVK4ddn7TRGeTSGktyGL0W9HOr6v1VdX1V3VhVH6mq05t1JpO8tmmpzyV5\nTZIDk9wR+Fazq6uTfCrJvwDHAR9qWrkHJjksydlJrkhyeZL/lSRLanh6kkuSXJXko0lut0ytH0ny\n3N2W7UjymGb6NUl2JrkmyQVJ7tLPj2DPh8pfJvlRkn9N8oilP6/lziXJf00yk2Q+yQ+SvGPJdt0k\nxzXTj0pycfPzubz5RUoaCwa5tPZOpNeift8K67wUuC9wPHCPZvqlVfUd4K7NOodX1UOq6g7A94BH\nN63cG4BzgAV6AX9P4KHAHwIkOQU4HXgMcGt6reRdAbibdwCnLs40QX074MNJHgbcH7hDVR0OPB64\nahU/h6V+k14Pwy2BvwTetOS7Zc8F+F/Ax6tqM9ABti/Zbml3/dnAM5tegbsB5+5jnVLrGOTS2rsl\n8MOq6q6wzqnAy6rqqqq6il4L/snNd9ntT5bOJzkSeCTwwqr6aVX9EHgt8IRmvWcD/6eq/rmp4ZXA\nCUl+dQ91vBe4x5LvTgXe0/yycANwKHCXJKmqb1fVzr5+Ar/s0qp6c/Ve7nAOcJskR/ZxLjcAxyY5\npqoWquqfdv95NBaAuyY5tKquqaod+1in1DoGubT2rgJulWSl/79uS6+Vvegy4DbN9N4Ghh0LHAj8\ne9NVfTXwN/Ra34vfv6757kdNPQUcs/uOquonwEf4RXD+PvD3zXefAf6a3rX9nUn+Jskhe6ltOVcu\nOeb1zeQhfZzLi+n9O/WVJBcledoy+//vwKOBy5J8Jsn99rFOqXUMcmntfRH4Gb2u7eV8n16ILToW\nuKLP/V8O/BS4ZVX9SlUdUVWbq+r45vvvAc9uvlv8/pCq+tIy+3sHcGoTfgc1AQ5AVf11Vd0buAtw\nJ3rBupZWPJeq+kFVPauqjgH+CDhr8br4UlV1XlUtXkp4P/DONa5TGlkGubTGqupaYCvwhiSnJLl5\nkk1JHplkcRT3PwAvTXKrJLcC/ifwtiW72dOAscX9Xwl8AnhNkkPTc1ySBzSr/C1wxuLAtCSHJ3nc\nCiV/hN4vEn8B/OOuApJ7J7lvczvc9fQCd6XLBau2t3NJ8rgkiz0J883xb1JDM/jv1CSHVdWNwI+B\nG9eyTmmUGeTSAFTVmcCL6A1q+wG9VvJz+cUAuJcDXwMuBC5opv/30l3svsvd5p8MTNK7Je1HwLvo\n3bZGVb2P3nXxf0gy3xzjESyjqhaA9wAPBt6+5KvDgDc2+/834If0BqqtpN/7xZeut+y5APcBvpzk\nWno/uxcsuXd86T6eBPxbc77PYskAPmmjS2/syYAP0rtWeB5weVWdnGQr8Ex6/8ABnFFVHxt4IZIk\nbTCb1uk4pwEX0/sNf9GZTatFkiTto4F3rSfpAI+id5/nTb4a9LElSdro1uMa+WvojXTdvQ//+c0T\npM5Ocvg61CFJ0oYz0K71JI8GdlbVjiTTS746C/iLqqokL6f3TOpn7GH7wV/AlyRphFTVqnqsB90i\nPwk4Ocl36d2r+ttJ3lpV/1G/GGX3RnojU/eoqsb2s3Xr1qHX4Ll7/p6/5+/5r99nXww0yKvqjKq6\nXVUdR+/JUedW1ZOTHL1ktccC3xhkHZIkbVTrNWp9d69OcgK9BztcSu/Z0JIkaZXWLcir6rPAZ5vp\nJ+9ldQHT09PDLmFoxvncwfP3/KeHXcJQjfv5r9a6PBBmX/VeuDS69UmStJaSUCM22E2SJA2QQS5J\nUosZ5JIktZhBLklSixnkkiS1mEEuSVKLGeSSJLWYQS5JUosZ5JIktZhBLklSixnkkiS1mEEuSVKL\nGeSSJLWYQS5JUosZ5JIktZhBLklSixnkkiS1mEEuSVKLGeSSJLWYQS5JUosZ5JIktZhBLklSixnk\nkiS1mEEuSVKLGeSSJLWYQS5JUosZ5JIktZhBLklSixnkkiS1mEEuSVKLGeSSJLWYQS5JUosZ5JIk\ntdimYRcgSdJK5ufnmZ2dZW5ujk6nw9TUFJs3bx52WSMjVTXsGpaVpEa5PknSYM3Pz7N9+3YWFhbo\ndrtMTEwwOTnJli1bNmSYJ6Gqsppt7FqXJI2s2dnZXSEO0O12WVhYYHZ2dsiVjQ6DXJI0subm5naF\n+KJut8vc3NyQKho9BrkkaWR1Oh0mJm4aVRMTE3Q6nSFVNHoMcknSyJqammJycnJXmC9eI5+amhpy\nZaPDwW6SpJE2TqPW92Ww27oEeZIJ4GvAXFWdnOQI4B+BY4FLgcdX1TV72M4glySNjVEetX4acMmS\n+dOBT1XVnYBzgZesUx2SJG0oAw/yJB3gUcDZSxafApzTTJ8DPGbQdUiStBGtR4v8NcCLgaV95EdV\n1U6AqroSOHId6pAkacMZ6CNakzwa2FlVO5JMr7DqshfCt23btmt6enqa6emVdiNJUnvMzMwwMzOz\nX/sY6GC3JK8Angj8HLg5cCjwXuDewHRV7UxyNPCZqvr1PWzvYDdJ0tgYucFuVXVGVd2uqo4DngCc\nW1VPAj4IPLVZ7SnA+wdZhyRJG9WwHgjzSuChSb4NPLiZlyRJq+QDYSRJGhEj17UuSZIGyyCXJKnF\nDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxy\nSZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrMIJckqcUMckmS\nWswglySpxQxySZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrM\nIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrMIJckqcUGGuRJDkry5STnJ7k4ySua5VuT\nzCX5evN5xCDrkCRpo0pVDfYAyS2q6rokBwBfAP4EeAjw46o6cy/b1qDrkyRpVCShqrKabQbetV5V\n1zWTBzXHu7qZX1WhkiTplw08yJNMJDkfuBKYqapLmq+en2RHkrOTHD7oOiRJ2ojWo0Xerap7Ah3g\nAUkeCJwFHFdVJ9AL+BW72CVJ0p5tWq8DVdW1ST4M3LuqPrvkqzcCH1xuu23btu2anp6eZnp6elAl\nSpK0rmZmZpiZmdmvfQx0sFuSWwE3VNU1SW4OfBx4GXBxVV3ZrPNC4D5VdeoetnewmyRpbOzLYLdB\nt8hvA5yTJPS68d9WVZ9O8tYkJwBd4FLg2QOuQ5KkDWngt5/tD1vkkqRxMpK3n0mSpMExyCVJajGD\nXJKkFjPIJUlqMYNckqQWM8glSWoxg1ySpBYzyCVJajGDXJKkFjPIJUlqMYNckqQWM8glSWoxg1yS\npBYzyCVJajGDXJKkFjPIJUlqMYNckqQWM8glSWoxg1ySpBYzyCVJajGDXJKkFjPIJUlqMYNckqQW\nM8glSWoxg1ySpBYzyCVJarFNwy5AkiSA+fl5ZmdnmZubo9PpMDU1xebNm4dd1shLVQ27hmUlqVGu\nT5K0Nubn59m+fTsLCwt0u10mJiaYnJxky5YtYxXmSaiqrGYbu9YlSUM3Ozu7K8QBut0uCwsLzM7O\nDrmy0WeQS5KGbm5ubleIL+p2u8zNzQ2povYwyCVJQ9fpdJiYuGkkTUxM0Ol0hlRRexjkkqShm5qa\nYnJycleYL14jn5qaGnJlo8/BbpKkkeCo9X0b7GaQS5I0Ihy1LknSmDHIJUlqMYNckqQWM8glSWox\ng1ySpBYzyCVJajGDXJKkFhtokCc5KMmXk5yf5OIkr2iWH5HkE0m+neTjSQ4fZB2SJG1UA38gTJJb\nVNV1SQ4AvgD8CXAycFVVvTrJnwFHVNXpe9jWB8JIksbGSD4QpqquayYPao53NXAKcE6z/BzgMYOu\nQ5KkjWjgQZ5kIsn5wJXATFVdAhxVVTsBqupK4MhB1yFJ0ka0adAHqKoucM8khwEfTzIN7N5fvmz/\n+bZt23ZNT09PMz09vfZFSpI0BDMzM8zMzOzXPtb1pSlJ/idwPfAMYLqqdiY5GvhMVf36Htb3Grkk\naWyM3DXyJLdaHJGe5ObAQ4HzgQ8AT21Wewrw/kHWIUnSRjXorvXbAOckCb1fGt5WVZ9urpm/M8nT\ngcuAxw+4DkmSNiTfRy5J0ogYua51SZI0WAa5JEktZpBLktRiBrkkSS1mkEuS1GIGuSRJLWaQS5LU\nYga5JEktZpBLktRiBrkkSS1mkEuS1GIGuSRJLWaQS5LUYga5JEkt1tf7yJMcABy1dP2q+t6gipIk\nSf3Za5An2QJsBXYC3WZxAccPsC5JktSHVNXKKyT/AvxmVV21PiXd5Ni1t/okSdooklBVWc02/Vwj\nvxy4Zt9KkiRJg9TPNfLvAjNJPgz8bHFhVZ05sKokSVJf+gny7zWfyeYjSZJGxF6vke9aMblFVV03\n4Hp2P6bXyCVJY2Mg18iTnJjkEuBbzfw9kpy1jzVKkqQ11M9gt9cCDweuAqiqC4AHDLIoSZLUn76e\n7FZVl++26MYB1CJJklapn8Fulyf5LaCSHAicBnxzsGVJkqR+9NMi/yPgecAxwPeBE4DnDrIoSZLU\nn35a5Heqqj9YuiDJScAXBlOSJEnqVz8t8u19LpMkSets2RZ5khOB3wJuneRFS746DDhg0IVJkqS9\nW6lrfRI4pFnn0CXLrwUeN8iiJElSf5YN8qr6LPDZJNdX1auXfpfk94DvDLo4SZK0sn6ukT9hD8te\nstaFSJKk1VvpGvkjgUcBxyR5/ZKvDgN+PujCJEnS3q10jfwK4GvAycB5S5b/GHjhIIuSJEn92evb\nz5JsqqqhtMB9+5kkaZzsy9vPVupaf2dVPR44P8kvpWlVHb8PNUqSpDW0bIs8yW2q6t+THLun76vq\nsoFWhi1ySdJ42ZcW+Ypd60kOAD5VVQ/a3+L2hUEuSRon+xLkK95+VlU3At0kh+9XZZIkaSD6eWnK\nT4CLknwS+M/FhVX1goFVJUmS+tJPkL+n+UiSpBHTT5DPAf9UVdevdudJOsBbgaOALvB3VbU9yVbg\nmcAPmlXPqKqPrXb/kiSNu37uIz8HOBH4ETALfA74fFVdvdedJ0cDR1fVjiSH0HuwzCnA/wB+XFVn\n7mV7B7tJksbGmt5HvqiqntLs/Lb03nr2BuC2fW57JXBlM/2TJN8EjlmsdzWFSpKkX7bXl6YkeWKS\nvwXeDTwE+GtgarUHSnJ74ATgy82i5yfZkeRsR8VLkrRv+rlG/lrgX4G/AT5TVZeu9iBNt/q7gdOa\nlvlZwF9UVSV5OXAm8Iw9bbtt27Zd09PT00xPT6/28JIkjaSZmRlmZmb2ax97vUYOkOSuwAOA+wN3\nBL5dVU/q6wDJJuBDwEer6nV7+P5Y4IN7euSr18glSeNkzR8I0+z0MOB2wLHA7YHD6Y1A79ebgUuW\nhngzCG7RY4FvrGJ/kiSp0c+o9QuBzzefz1XVXN87T06iN8r9IqCazxnAqfSul3eBS4FnV9XOPWxv\ni1ySNDbW/Fnrw2aQS5LGyUC61iVJ0ugyyCVJajGDXJKkFutn1PqrkxyW5MAkn07yH0meuB7FSZKk\nlfXTIn9YVV0L/Dd6I8zvALx4kEVJkqT+9BPki09/ezTwrqq6ZoD1SJKkVejnEa0fSvIt4HrgOUlu\nDfx0sGVJkqR+9PuI1l8BrqmqG5PcAjisebPZYIvzPnJJ0hhZ09eYJvntqjo3yWOXHmCJ96y+REmS\ntJZW6lp/IHAu8Dt7+K4wyCVJGjof0SpJ0ogY1NvPTmvuI0+Ss5N8PcnD9r1MSZK0Vvq5/ezpzX3k\nDwNuCTwJeOVAq5IkSX3pJ8gXm/iPAt5aVRcvWSZJkoaonyA/L8kn6AX5x5McSu894pIkacj2Otgt\nyQRwAvDdqppPckvgmKq6cODFOdhNkjRGBvU+8gLuArygmT8YuNkqa5MkSQPQT5CfBZwI/H4z/2Pg\nDQOrSJIk9a2fZ63/ZlXdK8n5AFV1dZLJAdclSZL60E+L/IYkB9DrYqd5aYqD3SRJGgH9BPnrgfcC\nRyb538DngVcMtCpJktSXft9+dmfgwfTuH/90VX1z0IU1x3XUuiRpbOzLqPV+g/wA4CiWXFOvqu+t\nusJVMsglSeNkTV9jumSnW4CtwE7gRnqt8gKO35ciJUnS2unngTD/Qm/k+lXrU9JNjm2LXJI0Ngb1\nQJjLgWv2rSRJkjRIy7bIk7yombwrcCfgw8DPFr+vqjMHXpwtcknSGFnra+SHNn9+r/lMNh9o7imX\nND7m5+eZnZ1lbm6OTqfD1NQUmzdvHnZZ0tjr5xr571XVu/a2bBBskUujYX5+nu3bt7OwsEC322Vi\nYoLJyUm2bNlimEtraFDXyF/S5zJJG9Ts7OyuEAfodrssLCwwOzs75MokLdu1nuSR9N5BfkyS1y/5\n6jDg54MuTNLomJub2xXii7rdLnNzc0OqSNKilVrkVwBfA34KnLfk8wHg4YMvTdKo6HQ6TEzc9J+L\niYkJOp3OkCqStGjZFnlVXQBckOTtVXXDOtYkacRMTU2xY8eOX7pGPjU1NezSpLHX1yNah8XBbtLo\ncNS6NHgDe9b6sBjk0vAY3NL6W9MgT/K2qnpSktOq6nVrUuEqGeTScHi7mTQca3372W8kuS3w9CRH\nJPmVpZ/9K1XSKPN2M6k9Vnqy298AnwaOozdafelvCNUsl7QBebuZ1B7Ltsir6vVV9evAm6vquKr6\nL0s+hri0gXm7mdQee32yW1U9J8k9kjy/+fgecmmDm5qaYnJycleYe7uZNLr6edb6C4BnAe9pFv0u\n8HdVtX2vO086wFuBo4Au8Maqen2SI4B/BI4FLgUeX1W/9KpUB7tJw+OodWn9DeT2syQXAidW1X82\n8wcDX6yqvbbMkxwNHF1VO5IcQu9a+ynA04CrqurVSf4MOKKqTt/D9ga5JGlsDOqlKQFuXDJ/Izcd\n+LasqrqyqnY00z8Bvgl06IX5Oc1q5wCP6bdgSZL0CyuNWl/0FuDLSd7bzD8GeNNqD5Tk9sAJwJeA\no6pqJ/TCPsmRq92fJEnqI8ir6swkM8D9m0VPq6rzV3OQplv93cBpVfWTJLv3ly/bf75t27Zd09PT\n00xPT6/m0JIkjayZmRlmZmb2ax8Df0Rrkk3Ah4CPLj4hLsk3gemq2tlcR/9Mc6vb7tt6jVySNDYG\ndY18f70ZuGS3x7x+AHhqM/0U4P3rUIckSRvOQFvkSU4CPgdcRK/7vIAzgK8A7wR+FbiM3u1n83vY\n3ha5JGlsDOr2s4OB66uqm+TXgDvT6yYf+DvKDXJJ0jgZVNf654CbJTkG+ATwJOD/WX15kiRprfV1\nH3lVXQc8Fjirqn4PuOtgy5IkSf3oK8iTnAj8AfDhZtkBgytJkiT1q58g/2PgJcB7q+riJMcBnxls\nWZIkqR99j1pPcoumi33dONhNkjROBjLYLcmJSS4BvtXM3yPJWftYoyRJWkP9dK2/Fng4cBVAVV0A\nPGCQRUmSpP709WS3qrp8t0U37nFFSZK0rvp5+9nlSX4LqCQHAqfRex2pJEkasn5a5H8EPA84Bvg+\nvVeRPm+QRUmSpP4M/O1n+8NR65KkcbIvo9b32rWe5PV7WHwN8LWq8q1lkiQNUT9d6zej153+neZz\nPNABnpHktQOsTZIk7UU/bz/7EnBSVd3YzG8CZoH7AxdV1V0GVpxd65KkMTKot58dARyyZP5g4Fea\nYP/Zag4mSZLWVj+3n70a2JFkBgi9h8G8onlP+acGWJskSdqLvkatJ7kNcN9m9qtVdcVAq/rFce1a\nlySNjX3pWu83yI8A7khv4BsAVfW5VVe4Sga5JGmcDOr2sz+k9zS3DrADuB/wReC396VISZK0dvoZ\n7HYacB/gsqp6EHBPYH6gVUmSpL70E+Q/raqfAiQ5qKq+BdxpsGVJkqR+9DNqfS7JZuB9wCeTXA1c\nNtiyJElSP1b1rPUkDwQOBz5WVQsDq+oXx3OwmyRpbKzpqPUkN6P35rM7ABcBb6qqn+93latgkEuS\nxslaP9ntHODe9EL8kcBf7UdtkiRpAFZqkV9UVXdvpjcBX6mqe61rcbbIJUljZK1b5DcsTqx3l7ok\nSerPSi3yG4H/XJwFbg5c10xXVR028OJskUuSxsiaPtmtqg7Y/5IkSdIg9fNAGEmSNKIMckmSWswg\nlySpxQxySZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklps2Zem\nrIUkbwL+G7Czqo5vlm0Fngn8oFntjKr62CDrkNRu8/PzzM7OMjc3R6fTYWpqis2bNw+7LGkkLPsa\n0zXZeXJ/4CfAW3cL8h9X1Zl9bO9rTKUxNz8/z/bt21lYWKDb7TIxMcHk5CRbtmwxzLXh7MtrTAfa\ntV5Vnweu3sNXqypS0vianZ3dFeIA3W6XhYUFZmdnh1yZNBqGdY38+Ul2JDk7yeFDqkFSC8zNze0K\n8UXdbpe5ubkhVSSNloFeI1/GWcBfVFUleTlwJvCM5Vbetm3brunp6Wmmp6cHXZ+kEdLpdLjiiitu\nEuYTExN0Op0hViWtjZmZGWZmZvZrHwO9Rg6Q5Fjgg4vXyPv9rvnea+TSmPMaucbJvlwjX48gvz29\nsL57M390VV3ZTL8QuE9VnbrMtga5JEeta2yMXJAneTswDdwS2AlsBR4EnAB0gUuBZ1fVzmW2N8gl\nSWNj5IJ8fxnkkqRxMnK3n0mSpMEyyCVJajGDXJKkFjPIJUlqMYNckqQWM8glSWoxg1ySpBYzyCVJ\najGDXJKkFjPIJUlqMYNckqQWM8glSWoxg1ySpBYzyCVJajGDXJKkFjPIJUlqMYNckqQWM8glSWox\ng1ySpBYzyCVJajGDXJKkFts07AKkjWR+fp7Z2Vnm5ubodDpMTU2xefPmYZclaQNLVQ27hmUlqVGu\nT1pqfn6e7du3s7CwQLfbZWJigsnJSbZs2WKYS+pLEqoqq9nGrnVpjczOzu4KcYBut8vCwgKzs7ND\nrkzSRmaQS2tkbm5uV4gv6na7zM3NDakiSePAIJfWSKfTYWLipv9LTUxM0Ol0hlSRpHFgkEtrZGpq\nisnJyV1hvniNfGpqasiVSdrIHOwmrSFHrUvaH/sy2M0glyRpRDhqXZKkMWOQS5LUYga5JEktZpBL\nktRiBrkkSS1mkEuS1GIGuSRJLWaQS5LUYga5JEktZpBLktRiBrkkSS020CBP8qYkO5NcuGTZEUk+\nkeTbST6e5PBB1iBJ0kY26Bb5W4CH77bsdOBTVXUn4FzgJQOuQZKkDWugQV5Vnweu3m3xKcA5zfQ5\nwGMGWYMkSRvZMK6RH1lVOwGq6krgyCHUIEnShjAKg9184bgkSfto0xCOuTPJUVW1M8nRwA9WWnnb\ntm27pqenp5menh5sdZIkrZOZmRlmZmb2ax+pGmyDOMntgQ9W1d2b+VcBP6qqVyX5M+CIqjp9mW1r\n0PVJkjQqklBVWdU2gwzKJG8HpoFbAjuBrcD7gHcBvwpcBjy+quaX2d4glySNjZEL8v1lkEuSxsm+\nBPkoDHaTJEn7yCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrMIJck\nqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnF\nDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxy\nSZJazCCXJKnFDHJJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFNg27ALXH/Pw8s7OzzM3N\n0el0mJqaYvPmzWNfiyQNU6pq2DUsK0mNcn3jZH5+nu3bt7OwsEC322ViYoLJyUm2bNmy7gE6SrVI\n0lpKQlVlNdvYta6+zM7O7gpOgG63y8LCArOzs2NdiyQN29C61pNcClwDdIEbquq+w6pFezc3N7cr\nOBd1u13m5ubGuhZJGrZhtsi7wHRV3dMQH32dToeJiZv+5zIxMUGn0xnrWiRp2IYZ5Bny8bUKU1NT\nTE5O7grQxevSU1NTY12LJA3b0Aa7JfkuMA/cCPxdVb1xD+s42G2EjNJI8VGqRZLWyr4MdhtmkN+m\nqv49ya2BTwLPr6rP77ZObd26ddf89PQ009PT61uoJEkDMjMzw8zMzK75l73sZe0J8psUkWwFflxV\nZ+623Ba5JGlstOb2syS3SHJIM30w8DDgG8OoRZKkNhvW7WdHAe9NUk0Nf19VnxhSLZIktdZIdK0v\nx651SdI4aU3XuiRJWhsGuSRJLWaQS5LUYga5JEktZpBLktRiBrkkSS1mkEuS1GIGuSRJLWaQS5LU\nYga5JEktZpBLktRiBrkkSS1mkEuS1GLDeo2pVjA/P8/s7Cxzc3N0Oh2mpqbYvHnzsMuSJI0gX2M6\nYubn59m+fTsLCwt0u10mJiaYnJxky5YthrkkbXC+xnQDmJ2d3RXiAN1ul4WFBWZnZ4dcmSRpFBnk\nI2Zubm5Gg3qzAAAJQ0lEQVRXiC/qdrvMzc0NqSJJ0igzyEdMp9NhYuKmfy0TExN0Op0hVSRJGmUG\n+YiZmppicnJyV5gvXiOfmpoacmWSpFHkYLcR5Kh1SRpP+zLYzSCXJGlEOGpdkqQxY5BLktRiBrkk\nSS1mkEuS1GIGuSRJLWaQS5LUYga5JEktZpBLktRiBrkkSS1mkEuS1GIGuSRJLWaQS5LUYga5JEkt\nZpBLktRiBrkkSS1mkEuS1GIGuSRJLWaQS5LUYga5JEktZpBLktRiQwvyJI9I8q0k/5zkz4ZVxyib\nmZkZdglDM87nDp6/5z8z7BKGatzPf7WGEuRJJoC/Bh4O3BX4/SR3HkYto2yc/2Me53MHz9/znxl2\nCUM17ue/WsNqkd8X+E5VXVZVNwD/AJwypFokSWqtYQX5McDlS+bnmmWSJGkVUlXrf9DkvwMPr6pn\nNfNPBO5bVS/Ybb31L06SpCGqqqxm/U2DKmQvvg/cbsl8p1l2E6s9GUmSxs2wuta/CtwhybFJJoEn\nAB8YUi2SJLXWUFrkVXVjkucDn6D3y8Sbquqbw6hFkqQ2G8o1ckmStDZG/sluSV6d5JtJdiT5/5Ic\nNuyaBm2cH5aTpJPk3CQXJ7koyQv2vtXGk2QiydeTjN0lpySHJ3lX8//9xUl+c9g1rZckL2nO+cIk\nf99cetywkrwpyc4kFy5ZdkSSTyT5dpKPJzl8mDUO0jLnv+rMG/kgp9f9fteqOgH4DvCSIdczUD4s\nh58DL6qquwInAs8bs/NfdBpwybCLGJLXAR+pql8H7gGMxWW3JMcCzwTuWVXH07v0+YThVjVwb6H3\nb91SpwOfqqo7Aeeysf/N39P5rzrzRj7Iq+pTVdVtZr9Eb4T7RjbWD8upqiurakcz/RN6/4iP1TMG\nknSARwFnD7uW9da0Pqaq6i0AVfXzqrp2yGWtl2uBBeDgJJuAWwBXDLekwaqqzwNX77b4FOCcZvoc\n4DHrWtQ62tP570vmjXyQ7+bpwEeHXcSA+bCcRpLbAycAXx5uJevuNcCLgXEcwPJfgB8meUtzaeHv\nktx82EWth6q6Gvgr4Hv0bsedr6pPDbeqoTiyqnZC7xd74Mgh1zNMfWXeSAR5kk8214QWPxc1f/7O\nknX+HLihqt4+xFK1TpIcArwbOK1pmY+FJI8Gdja9Emk+42QTcC/gDVV1L+A6el2tG16S44AXAscC\ntwUOSXLqcKsaCeP4C+2qMm9YD4S5iap66ErfJ3kqva7G316Xgoarr4flbGRNt+K7gbdV1fuHXc86\nOwk4OcmjgJsDhyZ5a1U9ech1rZc54PKq+loz/25gXAZ83hv4QlX9CCDJe4DfAsat8bIzyVFVtTPJ\n0cAPhl3Qeltt5o1Ei3wlSR5Br5vx5Kr62bDrWQc+LAfeDFxSVa8bdiHrrarOqKrbVdVx9P7uzx2j\nEKfpUr08ya81ix7M+Az6+zZwvyQ3SxJ65z4OA/1273n6APDUZvopwEb/Zf4m578vmTfy95En+Q4w\nCVzVLPpSVT13iCUNXPMX+Tp+8bCcVw65pHWT5CTgc8BF9LrUCjijqj421MKGIMkDgT+pqpOHXct6\nSnIPegP9DgS+Czytqq4ZblXrI8mL6YXYjcD5wB82g143pCRvB6aBWwI7ga3A+4B3Ab8KXAY8vqrm\nh1XjIC1z/mewyswb+SCXJEnLG/mudUmStDyDXJKkFjPIJUlqMYNckqQWM8glSWoxg1ySpBYzyKUB\nS3JUknck+U6Sryb5UJI77GWbv2weVfyqJLdK8qUk5zX32Ws3SbYmedGw65CGYSQe0SptcO8F3lJV\nvw+Q5O7AUcC/rLDNM4EjqqqSPAG4sKqeNfhSJbWNLXJpgJI8CFioqjcuLquqi6rqC833iy3vC5L8\nXrPs/cAhwHlJ/i/gVcBjmreBHZTkoUn+KcnXkvxjkls0290ryUzT6v9okqOa5S9IcnGSHc2TpHav\ncWJJHTuSPK9Z/uDmmBckOTvJgc3yf0uyrekhuGDxcapJHpDk/Gab85IcnOSBST645Fjbkzx5yX5e\n0Wzz1ab+jzc9F89ess2fJvlKU9vWJcv/PMm3k3wOuNPa/I1J7WOLXBqsuwHn7emLJI8Fjq+quyc5\nEvhqks9V1SlJrm3e/kWSncBvVNULktwSeCnw4Kq6vgn6FyV5JbCd3vOZr0ryeOAVwDPovXTk9lV1\nQ3rv+97ds+i9cev4pgdgc5KDgLcAD6qqf01yDvAc4PXNNj+oqt9I8hzgT5t9/Cnw3Kr6YvPLxU+b\ndVd6fOSlVXXPJGc2xzuR3nu4vwH8bZKHAnesqvs2zx//QJL703sr2uOB4+k9zvLrwNf2eARpgzPI\npeG5P/AOgKr6QZIZ4D7Ah1j+9aX3A+4CfKEJtgOBL9Jrkd4N+GSzfAK4otnmAuDtSd5H7znWu3sI\n8H9X87zmqppPcjzw3ar612adc4Dn8osgf2/z53nA7zbTXwBek+TvgfdU1fd7paxosbV+EXBwVV0H\nXJfkp80vHQ8DHprk683P5GDgjsBhwHubl0r8LMm4vVhI2sUglwbrYuBxfa67NPWWa8UG+ERV/cFN\nFiZ3A75RVXsaDPdo4AHAycCfJ7lbVXVXWc/uFt/KdCPNvyNV9aokH2qO94UkDwN+zk0v4d1smf10\nl0wvzm9qavg/Sy9NACQ5rY/6pbHgNXJpgKrqXGAyyR8uLkty96Z7eBb4H8016lsDU8CXF1dbZpdf\nAk5K8l+bfd0iyR3pvQLz1knu1yzflOQuzTa3q6rPAqfTa8kests+Pwk8O8kBzbZHNPs7NslxzTpP\nAmZWOtckx1XVxVX1anqv470zvbdX3SXJgUk203s1Zz8Wz//jwNOTHNwc47bNz+pz9MYNHJTkUOB3\n+tyvtOHYIpcG73eB1yU5HbgeuBT446r6fJIT6XV9d4EXV9V/NNvssUVeVT9M8lTgHc117AJeWlXf\nSfI4YHuSw4EDgNcm+Wfg/226qQO8rqqu3W23ZwO/BlyYZAF4Y1WdleRpwLubgP8q8Lcr1Qb8cXqD\n+26k1xPx0ea6/DvpXfP+N3rXsnedzgo/s8Vu/k8muTPwxaab/sfAE6vq/Ga/F9J7/eNXVtiXtKH5\nGlNJklrMrnVJklrMIJckqcUMckmSWswglySpxQxySZJazCCXJKnFDHJJklrs/we/N1SMapUmiQAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc635952c10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Defining the data\n",
    "coffee_v_thesis = np.array([\n",
    "        [0, 2.94],\n",
    "        [2.3, 5.45],\n",
    "        [1.2, 5.41],\n",
    "        [3.5, 12.10],\n",
    "        [5.1, 20.75],\n",
    "        [5.2, 21.48],\n",
    "        [10, 31.60],\n",
    "        [10.2, 40.84],\n",
    "        [6.2, 15.12],\n",
    "        ])\n",
    "\n",
    "# Visualising the data\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(x=coffee_v_thesis[:,0], y=coffee_v_thesis[:,1], color=\"#7a7a7a\", s=30)\n",
    "plt.ylabel(\"Pages of thesis written\")\n",
    "plt.xlabel(\"Coffees consumed\")\n",
    "plt.title(\"Coffee vs Thesis\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where P is the number of pages written and C is the number of coffees consumed: computing \n",
    "\n",
    "$$\\large{cov(P, C)}$$\n",
    "\n",
    "The documentation for `numpy.cov` can be found [here](https://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cov(P,C) = 43.66\n",
      "cov(P,C) (calculated using Numpy = 43.66\n",
      "**Notice the positive covariance, more coffee = more thesis pages!\n"
     ]
    }
   ],
   "source": [
    "\n",
    "C = coffee_v_thesis[:,0] # Select only the first column for coffee\n",
    "P = coffee_v_thesis[:,1] # Second column for pages\n",
    "n = coffee_v_thesis[:,0].shape[0] # Number of elements in the samples\n",
    "\n",
    "mean_C = np.mean(C)\n",
    "mean_P = np.mean(P)\n",
    "\n",
    "cov_P_C = np.sum((C - mean_C) * (P - mean_P)) / (n - 1)\n",
    "\n",
    "print (\"cov(P,C) = %0.2f\" % cov_P_C)\n",
    "\n",
    "# Computing the covariance using numpy\n",
    "# The first row, second element of the matrix corresponds to cov(P, C)\n",
    "cov_P_C = np.cov(coffee_v_thesis[:,0], coffee_v_thesis[:,1])[0][1]  \n",
    "\n",
    "print (\"cov(P,C) (calculated using Numpy = %0.2f\" % cov_P_C)\n",
    "\n",
    "print(\"**Notice the positive covariance, more coffee = more thesis pages!\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Covariance matrix\n",
    "$$\\large{\\overline{cov} =  \\begin{bmatrix}\n",
    "cov(X, X) & cov(X,Y) & cov(X,Z)\\\\ \n",
    "cov(Y, X) & cov(Y, Y) & cov(Y, Z)\\\\ \n",
    "cov(Z, X) & cov(Z, Y) & cov(Z, Z) \n",
    "\\end{bmatrix}}$$\n",
    "\n",
    "Computing the covariance matrix of `coffee_v_thesis`.  Notice the symmetry of the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The covariance matrix of coffee_v_thesis:\n",
      "[[  12.79027778   43.65606944]\n",
      " [  43.65606944  163.47413611]]\n"
     ]
    }
   ],
   "source": [
    "cov_P_C = np.cov(coffee_v_thesis[:,0], coffee_v_thesis[:,1])\n",
    "print(\"The covariance matrix of coffee_v_thesis:\\n%s\" % cov_P_C.__str__())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenvalues and Eigenvectors\n",
    "\n",
    "Calculate the eigenvalues and eigevectors of the covariance matrix `cov_P_C`\n",
    "\n",
    "There are a couple of different ways of calculating the eigenvalues and eigenvectors.  The most straight-forward is to use `numpy.linalg.eig`. The documentation for numpy.linalg.eig can be found [here](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigenvalues: [   1.05604092  175.20837297]\n",
      "eigenvectors: \n",
      "[[-0.96572297 -0.25957495]\n",
      " [ 0.25957495 -0.96572297]]\n"
     ]
    }
   ],
   "source": [
    "eigenvalues, eigenvectors = np.linalg.eig(cov_P_C)\n",
    "print(\"eigenvalues: %s\" % eigenvalues.__str__())\n",
    "print(\"eigenvectors: \\n%s\" % eigenvectors.__str__())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way of calculating the eigenvalues and eigenvectors is to use the singular value decomposition `numpy.linalg.svd` function within numpy (click [here](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.svd.html) for the documentation).  One convenient difference compared to `numpy.linalg.eig` is that singular value decomposition returnes the eigenvalues in ascending order which is useful to PCA\n",
    "\n",
    "The function `numpy.linalg.svd` takes a matrix `a` and returns two matrices `u` and `v` and an array of scalars `s` such that:\n",
    "\n",
    "`u * np.diag(s) * v = a` \n",
    "$$ a  = U S V^T $$\n",
    "\n",
    "Where `np.diag` computes the diagonal matrix of `s`.  \n",
    "\n",
    "The eigenvalues of the matrix `a` can be calculated as a multiple of `S`.\n",
    "\n",
    "The eigenvectors of the matrix `a` for the eigenvalues, correspond to the rows of `V` of the columns of `U`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('Correct transform: ', True)\n",
      "eigenvalues: [ 175.20837297    1.05604092]\n",
      "eigenvectors: \n",
      "[[-0.96572297 -0.25957495]\n",
      " [ 0.25957495 -0.96572297]]\n"
     ]
    }
   ],
   "source": [
    "U, S, V = np.linalg.svd(cov_P_C, full_matrices=False)\n",
    "\n",
    "eigenvalues = S\n",
    "eigevectors = U\n",
    "\n",
    "# You can check the decomposition occurred correctly if the following statement returns True\n",
    "print(\"Correct transform: \", np.allclose(cov_P_C, np.dot(U, np.dot(np.diag(S), V))))\n",
    "\n",
    "print(\"eigenvalues: %s\" % eigenvalues.__str__())\n",
    "print(\"eigenvectors: \\n%s\" % eigenvectors.__str__())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----------------------------------------------------------------------------------------------------------------------\n",
    "-----------------------------------------------------------------------------------------------------------------------\n",
    "# Executing PCA\n",
    "## 1. Subtract the means\n",
    "\n",
    "This ensures that the data is centered around the origin, looking at the graph below you can see the centre at (0,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Confirm the data is now centered about the mean: mu_x = -0.00, mu_y = 0.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7fc63585d910>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAYAAACmKP9/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcZHV97//XuxnajWUIKiodMcQVFZeruLaOa9QomkS9\nLlETkqsRM5p4zY1L7gNy7y9GzY3bGBNvNF40LlF/7ssVjUxocQFRFgGVqCANioo0aFAbqc/945we\na4bp7prpqqnq06/n41GPqTrn1KlPna6pd53v+Z7vSVUhSZK6ZWrcBUiSpOEz4CVJ6iADXpKkDjLg\nJUnqIANekqQOMuAlSeogA14bQpLnJvlekquTHJLkAUm+0T4+dtz17Ykkz0oyN4L1/jjJbfbieTdN\nckGSGwy7pnFK8u0kD23vvyTJ/56Amk5Jclx7/7FJ3j3umjS5DHitG0meluSMNoguTfKxJA8Y4Hmb\ngL8FHl5VB1XVlcBfAq9vH3941LXviSRfbX94XJ3kF0l+2r7nq5O8uF1s6ANYVNWBVXXRXjz1xcBb\nq+rnQy5pYlTVX1fVs9eyjiRHJOklGcr3blV9FDgqyV2GsT51jwGvdSHJC4FXA/8fcHPg1sDfAY8b\n4Om3AG4AXNA37Qjg/CGXORRVdZf2h8dBwBxwfBu+B1XVK8ZdX78k08CzgH8edy3rQGh+mGWI63w3\n8Jwhrk8dYsBr4iU5iGaP+/iq+lBV/bSqrquqj1fVi9tlppO8tt2zn0/ymiT7J7kd8LV2VVcm+XSS\nfweOBD7a7hXvn+SgJG9OclmSS5L8zyTpq+G4JOcnuSLJJ5LceplaP57k+F2mnZXkCe391yS5PMlV\nSc5OctQgm2D3L5W/SfKjJN9M8qj+7bXce0ny60m2J1lI8v0k7+p7Xi/Jke39xyQ5r90+l7Q/sHbn\nPsCVVXVZ33pOaV/ztLbl4UNJDk3yz+37/mL/9ktyxyQnt9v2giRP6pv3mCRfbp93cZIT+uYt7RE/\ns533/SQvXXYjrrCudv4zklyU5Ae7rifJCUne1t5/cJJLdpnf35x/77al6aok303yv9rF/q39d6Hd\nrvdpl1/2s5XkEe02uTLJNq7/WdgO/OZy71kbXFV58zbRN+A3gEVgaoVl/gfwOeDQ9nYa8JftvCOA\n64D0Lf9t4CF9jz8AvBG4IXBT4AvAf2nnPR74BnB7mh/FLwVOW6aOZwCf7Xt8FPAjYH/gkcAZwIHt\nvDsAh63y3k8Bjttl2rPa7XEczRf+HwGXDvhe3gm8pL0/Ddy/73nXAUe29y9bmgccDNx9mfqOBz6y\nm5q/AdwGOBA4r338kHb7nQS8pV32xsB3gGe27+VuwPeBO7bzHwTcub1/F+C7wLF9f9ce8Kb2vRwN\n/Ay4wzK1rrSuo4AfAw9o/1Z/227jh7bzTwDe1t5/MPCdXdb97b5lPwc8ve/9HbPC53DZz1b7t7sa\n+C1gP+BPgGv7Pw/AIe06Dxj3/1Nvk3dzD17rwaHAD6uqt8IyT6MJ9Cuq6gqaPf5ntvOyy7/0P05y\nc+DRwJ9W1c+q6ofAa4GntMs9B/jrqvpGW8MrgLsn+dXd1PEB4G59854GvL+qrqX5cj6Q5rhpqurr\nVXX5QFvg+i6qqn+qqqIJzFsmufkA7+Va4Igkh1fVYlV9btft0VoE7pzkwKq6qqrOWqaOzTTBuKu3\nVtVFVfVj4BPAhVV1Srv93gvco13uscC3q+pt1TgbeD/wJICqOrWqzmvvf5WmSfrBfa9TwIntezkH\nOJvmR8L1rLKu36H5oXJa+7f67+x9P4dF4LZJDq2qa6rq9F3m92/nlT5bjwa+WlUfqKbF6rXA93ZZ\n14/b9W3ey1rVYQa81oMrgJtm5c5Jt6LZE1xyMXDL9v5qX9RH0Oy1fbdt8r4S+AfgZn3zX9fO+1Fb\nTwGH77qiqvoJ8HF+GahPBd7RzjsFeANN34HLk/xDkgNWqW05O77oq+qn7d0DBngvf0bz//70JOcm\n+f1l1v87NE2/F7dN7vddZrkraX607Kr/h8tPd/N46X0fAdx3adu29T4NOAwgyX2SfKZtfl+gCcSb\nrvBa1/SteydJjllhXbcCdjS7V9U1NH/nvfEHNK0zX2sPR6zUhL7SZ2unmlq7Pj6wXX5hL2tVhxnw\nWg8+D/wceMIKy1xK82W55AiaZuZBXELTtHtoVf1KVR1SVZur6uh2/neA57TzluYfUFVfWGZ97wKe\n1obiDdpgB6Cq3lBV96JpEr4DTeAO04rvpaq+X1XPrqrDaZr237h03L1fVZ1ZVU+g+WHwIeA9y7ze\nOTTNy2upd/su2/agqvrjdv47gA8Ch1fVZprm+L3tpPbOFdb1XWBHi0ySG9O0HO3Of9A0vS8tux+/\n/AFFVX2zqp5WVTcDXgW8L8mN2P0PzZU+W9+l6Uzab9dWozvRtOb8ZKU3ro3JgNfEq6qraY6B/l2S\nxye5UZJNSR6dZKlX+buBv0hzTvZNaZpY3963mmVDoaq+B5wMvCbJgWkcmeRB7SJvAl6atkNckoOT\nPHGFkj9O8wPjfwD/sqOA5F7tXuQmmr3Yn9EcQx6a1d5LkicmWWp5WGhff6ca0nQ6fFqSg6rqOppm\n4OuWecnTgc1JbrnM/NV8FLh9kt9t/6b7t9vpDu38A2g68V2b5Biavfudyt2D11ppXe8DHpvk/kn2\np/nbLbfubwA3bD9/m4C/oOkD0BSUPL39DAJcRRPsPeAH7b+/3reulT5bH6M5nPOEJPsleQFty0af\nB9McApGux4DXulBVrwZeSPNl+n2aPZ/jafbIoDl97ks0e5Rnt/f/qn8Vu65yl8fPpPmSPp+mU9x7\naU6vo6o+SHNs9N1t0+45wKNYRlUt0hxHfhjNXuOSg4B/bNf/beCHwN+s+MYHPw7cv9yy7wW4N/DF\nJFfTbLvn1y/Pfe9fxzOAb7fv99lcP1ibJzTHq/9Pu/ye1rx0SOORNIc0Lmtvr6A5rRGav/H/THIV\nzd/+X3ZdxSqP+y27rqo6H3geTevLZTRN5fPL1Hx1u663tMv8eJdlHwWc127j1wD/uap+3h5K+Svg\ntLZJ/piVPlttX5InAa+k+az8Ok3n0X5PpfmRIF1Pmj46Y3rxZuSrU2m+jKaBD1XVS5McQvOf7wjg\nIuDJVXXV2AqVtKx2b/VU4B7V0cFukvwlTdP+H467liVJHgv8blU9ZdWFtSGNNeChOdZVVde0x7FO\nA/4rcCxwRVW9KsmfA4dUe76zJO1LSULTB+Gsqvqr1ZaXJsXYm+jb3qrQNMlN0fTKfTzNqT+0/67U\nuUqSRulMml7t/zjuQqQ9sWncBbSnPp1Jc3zpH6rq/CSHLZ0fXFXfa8/tlaR9rqruOe4apL0x9oBv\nB3e4R5rhSD+ZZAsDdpxJMt7jC5IkjUFVrXoGydib6Je0PVM/DtyLZhCQpYEubkHTa3q553kb4HbC\nCSeMvYb1cHM7ua3cTm6nSb8NaqwB356zfHB7/0bAI4CvAB8Gfq9d7Fk0A21IkqQBjbuJ/pbASW0v\n1Sng7VX1r0m+ArwnyXE0Q44+eZxFSpK03ow14KvqXOB6HViq6kfAw/d9Rd21ZcuWcZewLridBue2\nGozbaTBup+Eb+3nwa9FckGv91i9J0p5KQq2nTnaSJGl4DHhJkjrIgJckqYMMeEmSOsiAlySpgwx4\nSZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmS\nOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrI\ngJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICX\nJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySp\ngwx4SZI6yICXJKmDNo27AEmS1ruFhQXm5uaYn59nZmaG2dlZNm/ePNaaUlVjLWAtktR6rl+StP4t\nLCywbds2FhcX6fV6TE1NMT09zdatW0cS8kmoqqy2nE30kiStwdzc3I5wB+j1eiwuLjI3NzfWugx4\nSZLWYH5+fke4L+n1eszPz4+pooYBL0nSGszMzDA1tXOcTk1NMTMzM6aK2hrG+uqSJK1zs7OzTE9P\n7wj5pWPws7OzY63LTnaSJK3RvuxFP2gnOwNekqR1xF70kiRtYAa8JEkdZMBLktRBBrwkSR1kwEuS\n1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRB\nBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR001oBPMpPkM0nOS3Jukue30w9J\ncnKSryf5ZJKDx1mnJEnrTapqfC+e3AK4RVWdleQA4Ezg8cDvA1dU1auS/DlwSFW9eDfPr3HWL0nS\nvpaEqspqy411D76qvldVZ7X3fwJcAMzQhPxJ7WInAU8YT4WSJK1PE3MMPsltgLsDXwAOq6rLofkR\nANx8fJVJkrT+TETAt83z7wNe0O7J79rubju8JEl7YNO4C0iyiSbc315VH2onX57ksKq6vD1O//3l\nnn/iiSfuuL9lyxa2bNkywmolSdq3tm/fzvbt2/f4eWPtZAeQ5G3AD6vqhX3TXgn8qKpeaSc7SZJ+\nadBOduPuRf8A4FTgXJpm+AJeCpwOvAf4VeBi4MlVtbCb5xvwkqQNZV0E/FoZ8JKkjWZdnCYnSZJG\nw4CXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiA\nlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJck\nqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmD\nDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4\nSZI6yICXJKmDNo27AEmS1oOFhQXm5uaYn59nZmaG2dlZNm/ePO6ylpWqGncNey1Jref6JUnrw8LC\nAtu2bWNxcZFer8fU1BTT09Ns3bp1n4d8Eqoqqy1nE70kSauYm5vbEe4AvV6PxcVF5ubmxlzZ8gx4\nSZJWMT8/vyPcl/R6Pebn58dU0eoMeEmSVjEzM8PU1M6ROTU1xczMzJgqWp0BL0nSKmZnZ5ment4R\n8kvH4GdnZ8dc2fLsZCdJ0gAmpRf9oJ3sDHhJktYRe9FLkrSBGfCSJHWQAS9JUgcZ8JIkdZABL0lS\nBxnwkiR1kAEvSVIHGfCSJHWQAS9JUgcZ8JIkdZABL0lSBxnwkiR1kAEvSVIHGfCSJHWQAS9JUgcZ\n8JIkdZABL0lSBxnwkiR1kAEvSVIHGfCSJHWQAS9JUgcZ8JIkddCmQRZKsh9wWP/yVfWdURUlSZLW\nZtWAT7IVOAG4HOi1kws4eoR1SZKkNUhVrbxA8u/Afarqin1T0uCS1Gr1S5LUJUmoqqy23CDH4C8B\nrlp7SZIkaV8Z5Bj8t4DtST4G/HxpYlW9emRVSZKkNRkk4L/T3qbbmyRJmnCrHoPfsWBy46q6ZsT1\n7BGPwUuSNpqhHYNPcr8k5wNfax/fLckbh1Dj0vrfkuTyJOf0TTskyclJvp7kk0kOHtbrSZK0EQzS\nye61wG8AVwBU1dnAg4ZYw1vb9fd7MfDpqroD8BngJUN8PUmSOm+gkeyq6pJdJl03rAKq6rPAlbtM\nfjxwUnv/JOAJw3o9SZI2gkE62V2S5P5AJdkfeAFwwWjL4uZVdTlAVX0vyc1H/HqSJHXKIAH/R8Dr\ngMOBS4GTgeNHWdRuLNuT7sQTT9xxf8uWLWzZsmUflCNJ0r6xfft2tm/fvsfPG2QkuwdU1WmrTVuL\nJEcAH6mqo9vHFwBbquryJLcATqmqO+3mefailyRtKMMcyW7bgNPWIu1tyYeB32vvPwv40JBfT5Kk\nTlu2iT7J/YD7AzdL8sK+WQcB+w2rgCTvBLYAhyb5Ds2FbV4BvDfJccDFwJOH9XqSJG0EKx2DnwYO\naJc5sG/61cATh1VAVT1tmVkPH9ZrSJK00Swb8FX1b8C/JflpVb2qf16SJwEXjro4SZK0dwY5Bv+U\n3Uxz4BlJkibYSsfgHw08Bjg8yev7Zh0E/GLUhUmSpL230jH4y4AvAccCZ/ZN/zHwp6MsSpIkrc0g\n58FvqqqJ3GP3PHhJ0kYz6HnwKzXRv6eqngx8Jcn1UnRpUBpJkjR5lt2DT3LLqvpuO8rc9VTVxSOt\nbADuwUuSNppB9+BXbKJPsh/NZVsfMszihsWAlyRtNEMZqraqrgN6SQ4eWmWSJGnkBrma3E+Ac5N8\nCviPpYlV9fyRVSVJktZkkIB/f3uTJEnrxCABPw98rqp+OupiJEnScAxyHvxJwP2AHwFzwKnAZ6vq\nytGXtzI72UmSNpqh9KLfZYW3ormK3IuAW1XVIHv/I2XAS5I2mjUPdNO3ot8FZoG7Aj8E3kCzJy9J\nkibUIE30PwS+CfwDcEpVXbQP6hqIe/CSpI1mqE30Se4MPAh4IHA74OtV9Yw1V7lGBrwkaaMZykA3\n7YoOAm4NHAHcBjgY6K21QEmSNDqDNNGfA3y2vZ1aVfP7orBBuAcvSdpoht6LfhIZ8JKkjWZoTfSS\nJGn9MeAlSeogA16SpA4apBf9q5IclGT/JP+a5Aft4DeSJGlCDbIH/8iquhp4LHARcFvgz0ZZlCRJ\nWptBAn5pONvfBN5bVVeNsB5JkjQEg1ww5qNJvgb8FHhukpsBPxttWZIkaS0GHar2V4Crquq6JDcG\nDqqq7428utXr8jx4SdKGsuarySV5aFV9Jslv96+0z/vXVqIkSRqVlZroHwx8BnjcbuYVBrwkSRPL\noWolSVpHhnk1uRe058EnyZuTfDnJI4dTpiRJGoVBTpM7rj0P/pHAocAzgFeMtCpJkrQmgwT8UjPA\nY4C3VdV5fdMkSdIEGiTgz0xyMk3AfzLJgUBvtGVJkqS1WLWTXZIp4O7At6pqIcmhwOFVdc6+KHAl\ndrKTJG00w7wefAFHAc9vH98EuOEaapMkSSM2SMC/Ebgf8NT28Y+BvxtZRZIkac0GGYv+PlV1zyRf\nAaiqK5NMj7guSZK0BoPswV+bZD+apnrai83YyU6SpAk2SMC/HvgAcPMkfwV8Fnj5SKuSJElrMujV\n5O4IPIzm/Pd/raoLRl3YIOxFL0naaAbtRT9owO8HHEbfMfuq+s6aKhwCA16StNGs+XKxfSvaCpwA\nXA5cR7MXX8DRay1SkiSNxiAD3fw7TU/6K/ZNSYNzD16afAsLC8zNzTE/P8/MzAyzs7Ns3rx53GVJ\n69bQmuiTnAI8oqp+MazihsWAlybbwsIC27ZtY3FxkV6vx9TUFNPT02zdutWQl/bSmpvok7ywvfst\nYHuSjwE/X5pfVa9ec5WSOm1ubm5HuAP0ej0WFxeZm5vjcY973Jirk7ptpWPwB7b/fqe9Tbc3aM+J\nl6SVzM/P7wj3Jb1ej/n5+TFVJG0cywZ8Vf0lQJInVdV7++cledKoC5O0/s3MzHDZZZftFPJTU1PM\nzMyMsSppYxhkoJuXDDhNknYyOzvL9PQ0U1PNV83SMfjZ2dkxVyZ130rH4B9Ncw34w5O8vm/WQcDE\ndbiTNDn6e87f6U53AuAHP/iBveilfWilY/CXAV8CjgXO7Jv+Y+BPR1mUpPVr157zl112mT3npTEY\n5DS5/avq2n1Uzx7xNDlp8nzkIx/h9NNPv95x92OOOcae89IQDHqa3KrH4Cc13CVNJnvOS5NhkE52\nkjSwmZmZHZ3qlthzXtr3lg34JG9v/33BvitH0npnz3lpMix7DD7J+cDDgU8AW2guMrNDVf1o1MWt\nxmPw0mRy/HlpdNY8Fn2S5wPPBY4ELmXngK+qOnIYha6FAS9J2miGebGZv6+q5w6tsiEy4CVJG83Q\nAr5d2d2ApQNop1bVOWusbygMeEnSRjO00+Tapvp3ADdvb+9IsnXtJUqSpFEZpIn+HOB+VfUf7eOb\nAJ+vqqP3QX0rcg9ekrTRDG0PnqZz3XV9j69jlx71kiRpsqw0Fv2StwJfTPKB9vETgLeMriRJkrRW\ng3ayuyfwwPbhXFV9ZaRVDcgmeknSRjPUXvSTyoCXJG00wzwGL0mS1hkDXpKkDhrkPPibJJlq798+\nybFJ9h99aZIkaW8Nsgd/KnDDJIcDJwPPAP7PKIuSJElrM9B58FV1DfDbwBur6knAnUdbliRJWouB\nAj7J/YCnAx9rp+03upIkSdJaDRLwfwK8BPhAVZ2X5EjglNGWJUmS1mLg8+CT3Lhtqp8YngcvSdpo\nhnk1ufslOR/4Wvv4bkneOIQaJUnSiAzSRP9a4DeAKwCq6mzgQaMsSpIkrc1AA91U1SW7TLputwtK\nkqSJMMjV5C5Jcn+g2gFuXgBcMNqyJEnSWgyyB/9HwPOAw4FLgbu3jyVJ0oTyanKSJK0jg/aiX7WJ\nPsnrdzP5KuBLVfWhvSluUEkeRdPJbwp4S1W9cpSvJ0lSVwzSRH9Dmmb5C9vb0cAM8AdJXjuqwtoL\n3LyBpgf/nYGnJrnjqF5PkqQuGaST3dHAA6rqOoAkfw/MAQ8Ezh1hbccAF1bVxe3rvht4PO35+JIk\naXmD7MEfAhzQ9/gmwK+0gf/zkVTVOBzoPz1vvp0mSZJWMcge/KuAs5JsB0IzyM3Lk9wE+PQIaxvI\niSeeuOP+li1b2LJly9hqkSRp2LZv38727dv3+HkD9aJPckuaJnOAM6rqsj1+pT2U5L7AiVX1qPbx\ni4Hq72hnL3pJ0kYztLHoWz8DvgtcCdw2yb4YqvaM9rWOSDINPAX48D54XUmS1r1BTpP7Q5rR62aA\ns4D7Ap8HHjrKwqrquiR/DJzML0+TcwQ9SZIGsGoTfZJzgXsDX6iqu7enqr28qn57XxS4EpvoJUkb\nzTCb6H9WVT9rV3qDqvoacIe1FihJkkZnkF7080k2Ax8EPpXkSuDi0ZYlSZLWYo/Gok/yYOBg4P9W\n1eLIqhq8HpvoJUkbyqBN9MsGfJIb0lxJ7rY0I9a9pap+MdQq18iAlyRtNMM4Bn8ScC+acH808LdD\nqk2SJI3YSnvw51bVXdv7m4DTq+qe+7K41bgHL0naaIaxB3/t0p1Ja5qXJEkrW2kP/jrgP5YeAjcC\nrmnvV1UdtE8qXIF78JKkjWbQPfhlT5Orqv2GW5IkSdpXBh2LXpIkrSMGvCRJHWTAS5LUQYMMVStJ\nnbOwsMDc3Bzz8/PMzMwwOzvL5s2bx12WNDR7NFTtpLEXvaS9sbCwwLZt21hcXKTX6zE1NcX09DRb\nt2415DXxhnk1OUnqlLm5uR3hDtDr9VhcXGRubm7MlUnDY8BL2nDm5+d3hPuSXq/H/Pz8mCqShs+A\nl7ThzMzMMDW189ff1NQUMzMzY6pIGj4DXtKGMzs7y/T09I6QXzoGPzs7O+bKpOGxk52kDcle9Fqv\n1nw9+PXAgJckbTT2opckaQMz4CVJ6iADXpKkDjLgJUnqIANekqQOMuAlSeogA16SpA4y4CVJ6iAD\nXpKkDjLgJUnqIANekqQOMuAlSeogA16SpA4y4CVJ6iADXpKkDjLgJUnqIANekqQOMuAlSeogA16S\npA7aNO4CpI1gYWGBubk55ufnmZmZYXZ2ls2bN4+7LEkdlqoadw17LUmt5/q1MSwsLLBt2zYWFxfp\n9XpMTU0xPT3N1q1bDXlJeywJVZXVlrOJXhqxubm5HeEO0Ov1WFxcZG5ubsyVSeoyA14asfn5+R3h\nvqTX6zE/Pz+miiRtBAa8NGIzMzNMTe38X21qaoqZmZkxVSRpIzDgpRGbnZ1lenp6R8gvHYOfnZ0d\nc2WSusxOdtI+YC96ScMyaCc7A16SpHXEXvSSJG1gBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkd\nZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTA\nS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS\n1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EFjC/gkT0zy1STXJbnnLvNekuTCJBck\neeS4apQkab3aNMbXPhf4LeBN/ROT3Al4MnAnYAb4dJLbVVXt+xIlSVqfxrYHX1Vfr6oLgewy6/HA\nu6vqF1V1EXAhcMy+rk+SpPVsEo/BHw5c0vf40naaJEka0Eib6JN8CjisfxJQwMuq6iOjfG1Jkjay\nkQZ8VT1iL552KfCrfY9n2mm7deKJJ+64v2XLFrZs2bIXLylJ0mTavn0727dv3+PnZdx915KcAryo\nqs5sHx8FvAO4D03T/KeA3XayS2LfO0nShpKEqtq1/9r1jPM0uSckuQS4L/DRJJ8AqKrzgfcA5wMf\nB443xSVJ2jNj34NfC/fgJUkbzcTvwUuSpNEx4CVJ6iADXpKkDjLgJUnqIANekqQOMuAlSeogA16S\npA4y4CVJ6iADXpKkDhrpxWa0Pi0sLDA3N8f8/DwzMzPMzs6yefPmcZe1w6TXJ0mTwKFqtZOFhQW2\nbdvG4uIivV6Pqakppqen2bp160SE6KTXJ0mj5lC12itzc3M7whOg1+uxuLjI3NzcmCtrTHp9kjQp\nDHjtZH5+fkd4Lun1eszPz4+pop1Nen2SNCkMeO1kZmaGqamdPxZTU1PMzMyMqaKdTXp9kjQpDHjt\nZHZ2lunp6R0hunSMe3Z2dsyVNSa9PkmaFHay0/VMei/1Sa9PkkZp0E52BrwkSeuIveglSdrADHhJ\nkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6yICXJKmDDHhJkjrIgJckqYMMeEmSOsiAlySpgwx4SZI6\nyICXJKl7yZlWAAAJ/klEQVSDDHhJkjpo07gL6DqvXS5JGgevBz9CCwsLbNu2jcXFRXq9HlNTU0xP\nT7N161ZDXpK0V7we/ASYm5vbEe4AvV6PxcVF5ubmxlyZJKnrDPgRmp+f3xHuS3q9HvPz82OqSJK0\nURjwIzQzM8PU1M6beGpqipmZmTFVJEnaKAz4EZqdnWV6enpHyC8dg5+dnR1zZZKkrrOT3YjZi16S\nNEyDdrIz4CVJWkfsRS9J0gZmwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEG\nvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwk\nSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkd\nZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTAS5LUQQa8JEkdZMBLktRBBrwkSR1kwEuS1EEGvCRJHWTA\nS5LUQQa8JEkdZMBLktRBYwv4JK9KckGSs5L8/0kO6pv3kiQXtvMfOa4au2T79u3jLmFdcDsNzm01\nGLfTYNxOwzfOPfiTgTtX1d2BC4GXACQ5CngycCfg0cAbk2RsVXaE/3kG43YanNtqMG6nwbidhm9s\nAV9Vn66qXvvwC8BMe/9Y4N1V9Yuquogm/I8ZQ4mSJK1bk3IM/jjg4+39w4FL+uZd2k6TJEkDSlWN\nbuXJp4DD+icBBbysqj7SLvMy4J5V9Tvt423A56vqne3jNwMfr6r372b9oytekqQJVVWrHrreNOIC\nHrHS/CS/BzwGeGjf5EuBX+17PNNO2936PTYvSdJujLMX/aOAPwOOraqf9836MPCUJNNJfg24LXD6\nOGqUJGm9Guke/Cq2AdPAp9pO8l+oquOr6vwk7wHOB64Fjq9RHkeQJKmDRnoMXpIkjcek9KJfkyRb\n20Fxzk3yinHXM8mS/NckvSS/Mu5aJtVKgzCpObyW5GtJvpHkz8ddzyRKMpPkM0nOa7+Xnj/umiZd\nkqkkX07y4XHXMqmSHJzkve3303lJ7rPS8us+4JNsAR4H3LWq7gr8r/FWNLmSzACPAC4edy0TbreD\nMKn5EgbeAPwGcGfgqUnuON6qJtIvgBdW1Z2B+wHPczut6gU0h2a1vNfRnFV2J+BuwAUrLbzuAx54\nLvCKqvoFQFX9cMz1TLLX0HRs1ApWGIRJzaBTF1bVxVV1LfBu4PFjrmniVNX3quqs9v5PaL6IHc9j\nGe3Ox2OAN4+7lknVtiTOVtVbAdrB4K5e6TldCPjbAw9K8oUkpyS517gLmkRJjgUuqapzx13LOnMc\n8IlxFzFBdh2Iah6Da0VJbgPcHfjieCuZaEs7H3YKW96vAT9M8tb2UMb/TnKjlZ4wzl70A1thwJy/\noHkPh1TVfZPcG3gPcOS+r3L8VtlOL6Vpnu+ft2HtwSBM1y4NuiTtqSQHAO8DXtDuyWsXSX4TuLyq\nzmoPuW7o76YVbALuCTyvqr6U5LXAi4ETVnrCxFtpwJwkfwS8v13ujLYD2aFVdcU+K3BCLLedktwF\nuA1wdnvhnhngzCTHVNX392GJE2MvB2FSM+jUrfseLzsQ1UaXZBNNuL+9qj407nom2AOAY5M8BrgR\ncGCSt1XVM8dc16SZp2mF/VL7+H3Aip1cu9BE/0HaL+Ektwf234jhvpKq+mpV3aKqjqyqX6P5oNxj\no4b7alYYhElwBnDbJEckmQaeQjM4la7vn4Dzq+p14y5kklXVS6vq1lV1JM3n6TOG+/VV1eXAJW3O\nATyMVTolros9+FW8FfinJOcCPwf8YKyusBlsJbsdhGm8JU2GqrouyR/TnGkwBbylqlbsybsRJXkA\n8HTg3CRfofk/99Kq+r/jrUzr3POBdyTZH/gW8PsrLexAN5IkdVAXmuglSdIuDHhJkjrIgJckqYMM\neEmSOsiAlySpgwx4SZI6yICXxizJYUneleTCJGck+WiS267ynL9pL0P6yiQ3ba/FcGZ7/rV2keSE\nJC8cdx3SvtSFgW6k9e4DwFur6qkASe5KM07+v6/wnP9Ccw2GSvIU4JyqevboS5W0XrgHL41RkocA\ni1X1j0vTqurcqjqtnb+0p352kie10z4EHEBzPYH/BrwSeEJ7hakbJHlEks8l+VKSf0ly4/Z590yy\nvW0l+ESSw9rpz09yXpKzklzvwjpJpvrqOCvJ89rpD2tf8+wkb25H1yLJt5Oc2LYonL00tGaSByX5\nSvucM5PcJMmDk3yk77W2JXlm33pe3j7njLb+T7YtHc/pe86Lkpze1nZC3/SXJfl6klOBOwznLyat\nH+7BS+N1F+DM3c1I8tvA0VV11yQ3B85IcmpVPT7J1VV1z3a5y4H/VFXPT3IozdUDH1ZVP21/ALww\nyStohuA9tqquSPJk4OXAH9BcsOI2VXVtmmtO7+rZwBFtLZVkc5Ib0AwT/ZCq+maSk4DnAq9vn/P9\nqvpPSZ4LvKhdx4uA46vq8+2Pjp+1y640nOZFVXWPJK9uX+9+wI2BrwJvSvII4HZVdUx7IaUPJ3kg\ncA3wZOBommGHvwx8abevIHWUAS9NrgcC7wKoqu8n2Q7cG/goy19L4L7AUcBpbeDtD3yeZg/2LjTj\n64em9e6y9jlnA+9M8kGaizft6uHA31c7rnVVLSQ5GvhWVX2zXeYk4Hh+GfAfaP89E/it9v5pwGuS\nvAN4f1Vd2o71v5KlvftzgZtU1TXANUl+1v4YeSTwiCRfbrfJTYDbAQcBH2gvFvTzJF4QRxuOAS+N\n13nAEwdctj8Nl9vrDXByVT19p4nNJYO/WlW764T3m8CDgGOBlyW5S1X19rCeXS1dhe862u+Zqnpl\nko+2r3dakkcCv2DnQ4U3XGY9vb77S483tTX8df8hDoAkLxigfqnTPAYvjVFVfQaYTvKHS9OS3LVt\nZp4D/nN7DPxmwCzwxaXFllnlF4AHJPn1dl03TnI74OvAzZLct52+KclR7XNuXVX/BryYZs/3gF3W\n+SngOUn2a597SLu+I5Ic2S7zDGD7Su81yZFVdV5VvYrmsrN3BC4Gjkqyf5LNNJfAHMTS+/8kcFyS\nm7Svcat2W51K0y/hBkkOBB434HqlznAPXhq/3wJel+TFwE+Bi4A/qarPJrkfTRN6D/izqvpB+5zd\n7sFX1Q+T/B7wrvY4eQF/UVUXJnkisC3JwcB+wGuTfAP457a5O8DrqurqXVb7ZuD2wDlJFoF/rKo3\nJvl94H1t8J8BvGml2oA/SdOp8DqalotPtMf930NzTP3bNMfKd7ydFbbZ0uGCTyW5I/D5trn/x8Dv\nVtVX2vWeA1wOnL7CuqRO8nKxkiR1kE30kiR1kAEvSVIHGfCSJHWQAS9JUgcZ8JIkdZABL0lSBxnw\nkiR10P8DbeV8T+YyWx0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc6356745d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make a copy of the data\n",
    "X = np.copy(coffee_v_thesis)\n",
    "x_mean, y_mean = np.mean(X, axis=0)\n",
    "\n",
    "X[:,0] -= x_mean\n",
    "X[:,1] -= y_mean\n",
    "\n",
    "x_mean, y_mean = np.mean(X, axis=0)\n",
    "\n",
    "# Confirm the data is now centered about the mean i.e. the mean is close to 0\n",
    "print(\"Confirm the data is now centered about the mean: mu_x = %0.2f, mu_y = %0.2f\" % (x_mean, y_mean))\n",
    "\n",
    "# Visualising the data\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(x=X[:,0], y=X[:,1], color=\"#7a7a7a\", s=30)\n",
    "plt.ylabel(\"Pages of thesis written\")\n",
    "plt.xlabel(\"Coffees consumed\")\n",
    "plt.title(\"Coffee vs Thesis (mean adjusted)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Calculate the Covariance Matrix\n",
    "Calculate the covariance using the mean adjusted data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Covariance of mean adjusted data\n",
      "[[  12.79027778   43.65606944]\n",
      " [  43.65606944  163.47413611]]\n"
     ]
    }
   ],
   "source": [
    "cov_P_C = np.cov(X[:,0], X[:,1])\n",
    "print(\"Covariance of mean adjusted data\\n%s\" % cov_P_C.__str__())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Calculate the eigenvalues and eigenvectors\n",
    "Calculate the eigenvalues and eigenvectors of the covariance matrix.  We will use the singular value decomposition `np.linalg.svd` function to ensure the eigenvalue are ordered from largest to smallest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correct transform: True\n",
      "eigenvalues: [ 175.20837297    1.05604092]\n",
      "eigenvectors: \n",
      "[[-0.96572297 -0.25957495]\n",
      " [ 0.25957495 -0.96572297]]\n"
     ]
    }
   ],
   "source": [
    "U, S, V = np.linalg.svd(cov_P_C, full_matrices=False)\n",
    "\n",
    "eigenvalues = S\n",
    "eigevectors = U\n",
    "\n",
    "# You can check the decomposition occurred correctly if the following statement returns True\n",
    "print(\"Correct transform: %s\" % str(np.allclose(cov_P_C, np.dot(U, np.dot(np.diag(S), V)))))\n",
    "\n",
    "print(\"eigenvalues: %s\" % eigenvalues.__str__())\n",
    "print(\"eigenvectors: \\n%s\" % eigenvectors.__str__())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the percentage of the variance represented by each eigenvalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total variance V_t = 176.2644\n",
      "Proportion of variance V_prop = [ 0.99400877  0.00599123]\n",
      "Proportion of variance represented by principal component 99.40 %\n"
     ]
    }
   ],
   "source": [
    "# Calculate the total variance\n",
    "V_t = np.sum(eigenvalues)\n",
    "print(\"Total variance V_t = %0.4f\" % V_t)\n",
    "\n",
    "# Determine the proportion\n",
    "V_prop = eigenvalues / V_t\n",
    "print(\"Proportion of variance V_prop = %s\" % V_prop.__str__())\n",
    "\n",
    "# Determine the proportion represented by the principal component\n",
    "print(\"Proportion of variance represented by principal component %0.2f %%\" % (V_prop[0] * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will choose to use eigenvalues an eigenvectors that represent 98% of variance.\n",
    "\n",
    "$$\\large{P \\subset E \\ni\n",
    "P = E[..t]}$$\n",
    "\n",
    "In this example it is only the principal component that provides sufficient variance.  One option for PCA is to select only the principal component as it represents the most variance for any single eigenvalue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalue cumulative sum: [ 175.20837297  176.26441389]\n",
      "Cumulative sum  > f_v (0.99) = (array([0, 1]),)\n",
      "Number of selected eigenvalues = 1\n",
      "Selected eignvalues = [ 175.20837297]\n",
      "P = [[-0.96572297 -0.25957495]]\n"
     ]
    }
   ],
   "source": [
    "f_v = 0.99 # Proportion of total variation 98%\n",
    "\n",
    "#Select only the components that account for 98% of the variance\n",
    "cum_sum = np.cumsum(eigenvalues)\n",
    "selected_eigenvalues = np.nonzero(cum_sum > (f_v * V_t))\n",
    "print(\"Eigenvalue cumulative sum: %s\" % cum_sum.__str__())\n",
    "print(\"Cumulative sum  > f_v (%0.2f) = %s\" % (f_v, selected_eigenvalues))\n",
    "\n",
    "# Filter out the negligable components\n",
    "t = filter(lambda x: x!=0, selected_eigenvalues[0])[0]\n",
    "print(\"Number of selected eigenvalues = %d\" % t)\n",
    "print(\"Selected eignvalues = %s\" % str(eigenvalues[0:t]))\n",
    "# Assign the selected eigenvectors to P\n",
    "P = eigenvectors[0:t]\n",
    "print(\"P = %s\" % P.__str__())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Transform the data\n",
    "Now that we have the selected eigenvectors we can transform the mean adjusted data.  \n",
    "\n",
    "$$\\large{x_{transformed} = Px^T}$$\n",
    "\n",
    "The transformed data is no longer in terms of number of coffees consumed and pages written, but rather in terms of the selected vector.  The data lies in this new space and cannot be thought of directly in terms of the original parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transformed data: \n",
      "[[  8.41632937]\n",
      " [  5.54363342]\n",
      " [  6.61631168]\n",
      " [  2.65859245]\n",
      " [ -1.13188761]\n",
      " [ -1.41794962]\n",
      " [ -8.68031834]\n",
      " [-11.27193546]\n",
      " [ -0.7327759 ]]\n"
     ]
    }
   ],
   "source": [
    "data_transformed = P.dot(X.T).transpose()\n",
    "\n",
    "# Return the transformed data\n",
    "print(\"Transformed data: \\n%s\" % data_transformed.__str__())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice now the data set is one dimensional.  This is because we only selected one of the eigenvectors and the other axes has been \"thrown out\".  In this example there are only two eigenvectors, so by selecting both we would have included all the variance of the data set and not have gained much benefit from completing PCA.  We now have a representation of the original data set that is in terms of the *P*, the mean and is now smaller in size than the original."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Recovering the Data\n",
    "Now say we have done some processing on the data and we wish to return it to the original *coffees consumed* vs *pages written space*.  To recover the data we use:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAJoCAYAAADFzY2AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Wd4XNW59vH7GUkWLiO5N8ByqAnFASwChgMoFBMLTEli\nioPTg4HQSQEOBAEJgRwgvIeEA4QSG0IPwRRRA6ITIwVMMRgwccEG4YJtucmS5nk/7C0hyRppVEZb\nI/9/16VLs/szI2HdrLX22ubuAgAAQPeLRV0AAADAloogBgAAEBGCGAAAQEQIYgAAABEhiAEAAESE\nIAYAABARghiQAczsVDP7zMzWmNkgM9vfzD4Il4+Kur72MLMfmNmLaThvlZmN7cBxQ83sPTPL7eqa\nomRm/zGzg8PXF5jZzT2gpufM7Mfh6yPN7J6oawKiRhADuomZTTWz18PAsMTMHjOz/VM4LlvSNZIO\ndfc8d/9C0qWS/jdcfjjdtbeHmb0TBsQ1ZlZrZhvC97zGzM4Pd+vyCQzdPe7uCzpw6PmSbnf36i4u\nqcdw99+7+8mdOYeZFZhZwsy65O+Guz8qaRcz260rzgdkKoIY0A3M7FxJ10r6raThksZI+rOkySkc\nPlJSrqT3Gq0rkDS3i8vsEu6+WxgQ8yS9KOm0MCTlufuVUdfXmJn1kfQDSXdGXUsGMAUB2rrwnPdI\nmt6F5wMyDkEMSDMzy1PQgnWau89y9w3uXufupe5+frhPHzO7Lmwp+8TM/mhmOWa2o6T3w1N9YWbP\nmNlHkraT9GjYypRjZnlmdouZLTWzxWZ2uZlZoxp+bGZzzWyFmT1uZmOS1FpqZqc1W/emmR0Tvv6j\nmVWa2Wozm2Nmu6TyEbR8KfsfM1tpZvPN7FuNP69k78XMtjezMjNbZWafm9ndjY5LmNl24etiM3s3\n/HwWh0G4JftI+sLdlzY6z3PhNV8OW/JmmdkQM7szfN//avz5mdlXzeyp8LN9z8ymNNpWbGb/Do9b\naGaXNNpW38L0/XDb52Z2YdIPsZVzhdunmdkCM1vW/DxmdomZzQxfH2Rmi5ttb9yNuXfYcrvazD41\ns6vD3Z4Pv68KP9d9wv2T/m6Z2WHhZ/KFmV2vzX8XyiQdkew9A1sCghiQfhMUtGg91Mo+F0n6hqRx\nkr4evr7I3T+UtGu4T767H+ruO0haJOmIsJWpRtIMSZsUBLQ9JR0m6aeSZGZHK+h+O0bSMAWtVA0B\nppm7JU2tXwiD1hhJj5nZREn/JWkHd8+XdJykFe34HBrbR0EL3xBJ/yPp1kbbkr4XSZdLetLdB0ra\nRtL1jY5r3N15i6Sfha1yu0l6Nkkdu0ua18L64yV9T9JoSTtIejWscZCCYHyJJJlZP0lPKWhRGyrp\nBEl/NrOvhudZK2la+HkdIekU23xM3/6SdpR0qKTfmNnOSWpNeq7w53RDo5qHSNo6yXmk1ruG/5+k\n68LrbC/pvnD9geH3vPD37l+t/W6Z2VBJf5d0oYLPZn74Xht7T1KBmQ1opR6gVyOIAek3RNJyd0+0\nss9USZe6+wp3X6GgBe374TZr9l2Nl81suKRJks5x943uvlzSdQpCgRR0/fze3T8Ia7hS0h5mtm0L\ndfxD0tcbbZsq6cEw7NVIiisY12PuPs/dK1P6BDa3wN1v8+BhtzMkjTKz4Sm8lxoFf7i3dvdN7v5K\n888jtEnSrmYWd/fV7v5mkjoGSqpqYf3t7r7A3askPS7pQ3d/Lvz87lcQECXpSEn/cfeZHpgj6UFJ\nUyTJ3V9w93fD1+8o6Io7qNF1XFJJ+F7ekjRHQRDfTBvn+o6kR9z95fBndbE6Pg5vk6QdzGyIu693\n99nNtjf+nFv73Zok6R13/0fYAnydpM+anasqPN/ADtYKZDyCGJB+KyQNtdYHOY9W0MpVb6GkUeHr\ntv6gFkjKkfRp2NX3haQbFbRQ1G//f+G2lWE9rhZaTNx9raRSfRl8TpT0t3Dbc5L+pGBsW6WZ3diJ\nloyGP8juviF8OSCF9/JLBf9uzTazt83sR0nO/x0FrUYLw67GfZPs94WCcNlc44C5oYXl+vddIGnf\n+s82rHeqpBGSZGb7mNmzYbfjKgXBZWgr11rf6NxNmNk3WjnXaEkN3Y3uvl4db638iaSdJb0fdsO2\n1nXY2u9Wk5pCzZfj4f6rOlgrkPEIYkD6vSqpWkH3TTJLFPxRq1cgaWmSfZtbLGmjpCHuPtjdB7n7\nQHcfF25fJGl6uK1++wB3fy3J+e6WNDUML7lhAJMkufuf3L1Q0i4K/lj/MsUaU9Xqe3H3z939ZHff\nWtIpkm6oHxfWmLtXuHt9d9ksfdm91txbknbqZL1lzT7bPHc/Pdz+NwVd0luH3ak3qeOD3e9q5Vyf\nSmpo4Qy7TIckOc86Sf0a7ZulL4Ou3H2+u09192GS/iDpATPrq5b/h6C1361PFXRrN9a8FfZrClpH\n17b2xoHejCAGpJm7r1EwpujPZna0mfU1s2wzm2Rm9XcR3iPpIgvmtBqqoGvpjkanSfrH290/UzBO\n6Y9mFrfAdmZWP6bnJkkXhuOIZGb5ZvbdVkouVRAEL5N0b0MBZoVhq0y2glahjZJa625tt7bei5l9\n18zqW/JWhddvUoMFNy9MNbM8d69T0P1Vl+SSsyUNNLNRSba35VFJO5nZSeHPNCf8nOrHeQ1QcDNA\njZl9Q43G39WX245rtXauByQdaWb7mVmOgp9dsnN/IGmr8PcvW8H4xD4NBZl9L/wdlKTVCgJYQtKy\n8Pv2jc7V2u/WYwq6sY8xsywzO0thS2EjByno+gW2WAQxoBu4+7WSzlXwR+9zBS0Jp+nLAfy/lVSu\noIVmTvj6d41P0fyUzZa/r+CP6VxJKxWMYxoZXvshBWN37gm7tN6S9C0l4e6bFIxzOkRBK0y9PEl/\nCc//H0nLFQy0b02q45Qa75f0vUjaW9K/zGyNgs/uTP9y7rDG55gm6T/h+z1Zmweg4IBgPNVfw/3b\nW3N9V+5EBV25S8OvKxXcnCEFP+PLzWy1gp/9vc1P0cZyY0nP5e5zJf1cQWvmUgVdhJ8kqXlNeK5b\nw32qmu37LUnvhp/xHyUd7+7VYRfy7yS9HHZFfqO1361wrOMUSVcp+F3ZXtLLzco5UUGYA7ZYFoyV\nTfNFgrExFZIWu/tRFtx2/TMFf5Ak6UJ3fyLthQBAM2HrzwuS9vReOqmrmV2qoEvzp23u3E3M7EhJ\nJ7n7CW3uDPRi2d10nbMkvavg/6jrXRu2EgBAZMI7M1OZDy0jmZkpeH/J7hyNhAcz6z8adR1A1NLe\nNWlm20gqVjCvT5NN6b42AEAVCu5i/EvUhQDYXHeMEfujgjurmveBnm7BjN23mFl+N9QBAFscd9/L\n3fdz98/b3htAd0trEAvnn6kMJ1Ns3AJ2g6Tt3H0PBfMJ0UUJAAC2OGkdrG9mV0g6SVKtpL4KJu97\n0N2/32ifAgUzQo9r4fj030kAAADQRdy9XUOvuuWuSSl40Kyk88K7JkeG8wXJzM6RtLe7b3Z7efAU\nFbJYpiopKVFJSUnUZaAD+NllNn5+mYufXWYzs3YHse66a7K5P5jZHgomB1yg4FEdAAAAW5RuC2Lu\n/ryk58PX329jdwAAgF6PmfWRNkVFRVGXgA7iZ5fZ+PllLn52W55uGyPWEYwRAwAAmSKTxogBANAu\nY8eO1cKFC6MuA1BBQYEWLFjQJeeiRQwAkBHC1oaoywCS/i52pEWMMWIAAAARIYgBAABEhCAGAAAQ\nEYIYAAAR+f3vf6+TTz65y/dtSywW08cff9wl50LnEMQAAL3esmXL9J///EcbNmxI2zX++te/aty4\ncerfv79Gjx6t0047TatXr271mAsuuEA333xzSudvz75tMUs+nryoqEh9+/ZVfn6+Bg4cqL333ltX\nXXWVNm3alPL5Oxv0XnrpJe2///4aOHCghg4dqgMOOEAVFRXdcu3uxvQVAICMVldXp88//1zZ2dka\nOnRok5CxevVq3XHHHVq+fLlisZjq6uo0YcIEHX744a2Gkfa65pprdPXVV2vmzJk6+OCDtWTJEp16\n6qk67LDD9Morryg7e/M/t3V1dcrKyuqyGtqjtbtPzUw33HCDfvSjH2nDhg16/fXXddZZZ+npp5/W\nM888k9L5O/PZVlVVafLkybrppps0ZcoUbdq0SS+++KJyc3PTfu0o0CIGAMhYc+fO1ZVXXqmbb75Z\nf/7zn/XHP/5Rn376qaQgbNx2222qrKxUTU2NqqurVVtbq9dee02zZ89uOMeGDRv0/PPP6/bbb9eD\nDz6oJUuWtKuGqqoqlZSU6E9/+pMOO+wwZWVlacyYMbrvvvu0YMEC3XnnnZKkSy+9VFOmTNG0adM0\ncOBAzZgxQ5deeqmmTZvWcK6ZM2dq7NixGjZsmH7729/qK1/5ip599tmG4+v3XbhwoWKxmGbOnKmC\nggINHz5cV1xxRcN5Xn/9de23334aNGiQtt56a51xxhmqra1N+T3VB7W+ffvqwAMP1MMPP6xXX31V\npaWlbZ7/oIMOkrtr3LhxysvL0/33369Vq1Zp8uTJGj58uIYMGaLJkydr6dKlLV77gw8+kJnpuOOO\nk5kpNzdXhx56qHbbbbeGfW677TbtsssuGjJkiCZNmqTFixcnvXZPRxADAGSkJUuW6L777tP69eu1\nadMm1dTUaMWKFbrlllu0YcMGLVq0SGvWrFEikWhyXE1NjZ5//nlJQYvZddddp2effVYfffSR/v3v\nf+svf/mLXn/99ZTreOWVV1RdXa1jjz22yfr+/furuLhYTz/9dMO6hx9+WMcdd5xWrVqlqVOnSvqy\nBWfu3Ln6+c9/rrvvvluffvqpVq9evVlYad7a8/LLL+vDDz/UM888o8suu0zz5s2TJGVlZem6667T\nypUr9eqrr+rZZ5/VDTfckPJ7am7bbbdVYWGhXnzxxTbPX//Zvv3221qzZo2mTJmiRCKhH//4x1q8\neLEWLVqkfv366fTTT2/xWjvttJOysrL0wx/+UE888YRWrVrVZPusWbN05ZVX6qGHHtKyZct0wAEH\n6IQTTkh67Z6OIAYAyEgvvvhii608dXV1mjNnzmZ/wBtbt26dJOmJJ57QunXrGs7j7qqpqdGjjz6q\njRs3plTH8uXLNXToUMVim/9JHTVqlJYvX96wPGHCBE2ePFmStNVWWzXZ9+9//7uOOuooTZgwQdnZ\n2brssstava6ZqaSkRH369NG4ceP09a9/XXPmzJEk7bXXXvrGN74hM9OYMWN08sknN4SUjho9erRW\nrlyZ8vkbd38OHjxYxx57rHJzc9W/f39dcMEFSeuJx+N66aWXFIvFdPLJJ2v48OE6+uijtWzZMknS\nTTfdpAsuuEA77bSTYrGYzj//fL355psNrWLNr93TEcQAABlp2bJlLf7Bramp0fLlyzVq1Kikf5CH\nDh0qSXr//fdb3CcrK0vz589PqY6hQ4dq+fLlm7W8SdKnn37acC0paFlKZunSpU229+3bV0OGDGn1\n2iNGjGh43a9fP61du1aS9OGHH2ry5MkaNWqUBg4cqP/+7/9uEgg7YsmSJRo8eHCHzr9hwwZNnz5d\nY8eO1cCBA3XQQQdp1apVSX8+O++8s2677TYtWrRI77zzjpYuXaqzzz5bUtAte9ZZZ2nw4MEaPHiw\nhgwZIjNrd5dyT0EQAwBkpK233rrFgdl9+vTRyJEjNXz4cI0dO3azgfI5OTmaOHFim+dPddD3hAkT\nlJubqwcffLDJ+rVr1+rxxx/XoYcemtI5R40apU8++aRhecOGDVqxYkVKNTR36qmn6mtf+5rmz5+v\nVatW6Xe/+12nWokWL16siooKHXjggR06/zXXXKMPP/xQr7/+ulatWqUXXnhBUmotVzvttJN++MMf\n6p133pEUhNmbbrpJK1eu1MqVK/XFF19o7dq12nfffTv8/qJEEAMAZKQDDjhgs5BlZg1ddZL0ve99\nT3vttZeys7MVi8U0cOBATZkyRTvvvLMkadddd22xSzGRSGj77bdPqY68vDz95je/0RlnnKEnn3xS\ntbW1WrBggY4//niNGTNGJ510Ukrn+e53v6tHHnlEr732mmpqalRSUtLq/q2FmKqqKuXl5alfv356\n//339X//938p1dBc/Y0MxxxzjPbdd19NmjQppfOPHDmyyRQSVVVV6tu3r/Ly8rRy5cpW39u8efN0\n7bXXNrRwLV68WHfffbcmTJggSTrllFN0xRVXaO7cuZKCcX4PPPBA0mv3dAQxAEBGGjZsmH74wx9q\n2LBhysrKUlZWlgoKCjR9+nT16dNHUtD6dfTRR+s3v/mNLrroIv3iF7/Qrrvu2nCOb33rW4rH48rJ\nyZEUzEGVnZ2tY445JuXpEiTpl7/8pa644gr94he/UH5+viZMmKCCggI988wzDeduyy677KLrr79e\nxx9/vEaPHq28vDwNHz48aR3NW9caL1999dX629/+pry8PE2fPr1hMHuyY5s7/fTTlZ+fr5EjR+rc\nc8/VlClT9Pjjj6d8/pKSEn3/+9/X4MGD9cADD+icc87R+vXrNXToUO23334qLi5Oeu14PK5//etf\n2meffRSPx7Xffvtp3LhxuvrqqyVJxxxzjM4//3ydcMIJGjhwoMaNG6cnnngi6bV7OuvJA9rMzHty\nfQCA7mNmSVuB1q1bp6ysrM0GwKdi06ZNevPNNzV//nzl5+dr77331rBhwzpbbqetW7dOAwcO1Ecf\nfaSCgoKoy0EjyX4Xw/XtmsiMIAYAyAitBbHe4tFHH9UhhxyiRCKh8847T6+//nrKM8qj+3RlEKNr\nEgCAHmLWrFkaPXq0ttlmG82fP1/33HNP1CUhzWgRAwBkhC2hRQyZgRYxAACAXoAgBgAAEBGCGAAA\nQEQIYgAAABEhiAEAAESEIAYAANplxowZOuCAA6Iuo1cgiAEAer+qj6Tlr0mbVqfl9GPHjlW/fv2U\nl5enUaNGadq0aaqqqkrLtXqKZI9JWrhwoWKxmPLy8ho+j6OOOkrPPPNMyufuiqB3xRVXaLvttlNe\nXp7GjBmjE088sduu3R4EMQBAZkvUSqvnSlXzpeZzO234VHquWCqbLL32M+mJb0jvXLn5fp1kZnrs\nsce0Zs0azZkzR2+//bZ++9vfduk10qWurq7Lz2lmWr16dcPnceihh+rYY4/VzJkzUzre3dt8HmZr\nZsyYob/97W969tlntWbNGpWXl+uQQw7plmu3F0EMAJC5lj4pPbG39OIUqexI6ZlvBqFMCsLWy9+T\n1syT6jZItVVSolr6zwxpwZ1fnmPTaumDG6WXp0n//pW06p0OlVI/wefw4cN1+OGH69133/3yEps2\n6Re/+IUKCgo0atQonXbaaaqurm7YPmvWLO25557Kz8/XjjvuqKeeekqS9Omnn+roo4/WkCFDtNNO\nO+mWW25pWN+vXz+tWrWq4RxvvPGGhg0b1hCsbrvtNu2yyy4aMmSIJk2apEWLFjXsG4vFdMMNN2in\nnXbSTjvtJEl6//33NXHiRA0ZMkRf+9rXdP/99zfsv3LlSh111FHKz8/Xvvvuq/nz57fr8zjzzDNV\nUlKiX//61w3br7rqKu2www7Ky8vTbrvtpoceeqihjlNPPVWvvvqq4vG4Bg8eLEkqLS3VXnvtpfz8\nfBUUFOjSSy9Neu3y8nIdfvjhGjt2bEMNP/3pTxu2r1mzRj/96U81evRobbvttrr44ovl7kmvnU4E\nMQBAZlr1jlRxlrTpC6l2XRC21i2QXjw+CFcrK6QNlZI3a/Gp2yB98H/B6w2fSf88RHr/OmnZi9Ki\nB6QXvystuLvDZX3yySd6/PHHtc8++zSs+/Wvf62PPvpIb731lj766CMtWbJEl112mSRp9uzZ+sEP\nfqBrrrlGq1ev1gsvvNAQII4//niNGTNGn332me6//35deOGFKisr06hRo7Tffvvp73//e8M17r77\nbk2ZMkVZWVmaNWuWrrzySj300ENatmyZDjjggM265mbNmqXZs2dr7ty5Wr9+vSZOnKiTTjpJy5cv\n1z333KPTTjtN77//viTptNNOU79+/VRZWalbb71Vt912W7s/l29/+9v6/PPPNW/ePEnSDjvsoJdf\nfllr1qzRJZdcopNOOkmVlZX66le/qhtvvFETJkxQVVWVVq5cKUkaMGCA7rjjDq1evVqPPfaYbrzx\nRj388MMtXmvffffVzJkzdfXVV6uiokKJRKLJ9h/84Afq06ePPv74Y73xxht6+umndcsttyS9dlq5\ne4/9CsoDAMB9s78Js093/8dX3P9R0PTr4a+6z5/hvvgh90d22Xz7PwrcZ+0YnuMM939s18L2nd03\nrU65trFjx3o8Hvd4PO5m5sccc4zX1dU1bO/fv79//PHHDcuvvPKKf+UrX3F39+nTp/u555672TkX\nL17s2dnZvm7duoZ1F1xwgf/oRz9yd/dbbrnFDz744IZt2267rb/00kvu7j5p0iS/7bbbGrbV1dV5\nv379fNGiRe7ubmZeVlbWsP3ee+/1Aw88sMn1p0+f7pdddpnX1dV5Tk6Of/DBBw3bLrzwQj/ggANa\n/CwWLFjgsVisyft3d9+4caObmb/yyistHrfHHnv4ww8/7O7uf/3rX5Oev97ZZ5/d4udW76677vLD\nDjvMBwwY4EOHDvWrrrrK3d0rKys9NzfXN27c2LDv3Xff7d/85jdTvnayfBKub1fWoUUMAJCZquZL\nSmy+vm6DtPZjKW+XzVvD6g3YPvj+2TOSWtgnli0te6Vd5cyaNUtr1qxRWVmZnnvuOVVUVEiSli1b\npvXr12v8+PEaPHiwBg8erEmTJmnFihWSpMWLF2v77bff7HxLly7V4MGD1a9fv4Z1BQUFWrJkiSTp\nO9/5jl577TVVVlbq+eefV1ZWlvbff39JwYD5s846q+F6Q4YMkZk1HCtJ22yzTcPrhQsX6rXXXmvY\nf9CgQbrrrrtUWVmpZcuWqba2tsn+BQUF7fpsJDVcu767b+bMmdpzzz01aNAgDRo0SO+++66WL1+e\n9PjZs2fr4IMP1vDhwzVw4EDddNNNre5/4okn6qmnntKqVat044036uKLL9bTTz+thQsXqqamRqNG\njWp4r6ecckqr50qn7EiuCgBAZw3cLRj/1TxIZfWT8neR8naUhuwT3C2ZqG60fStp118Fr1sdlN2+\nAdsejok68MADdfrpp+tXv/qVnnvuOQ0dOlT9+vXTu+++q1GjRm123LbbbtvimKvRo0dr5cqVWrdu\nnfr37y9JWrRokbbeemtJ0sCBAzVx4kTdc889eu+993TCCSc0HDtmzBhddNFFrd4p2HhA+rbbbqui\noiI9+eSTm+2XSCSUk5OjxYsXN4wnazzeLFUPPvigRowYoZ133lmLFi3SySefrOeee04TJkyQJO25\n554Nn2FLg+WnTp2qM888U08++aRycnJ0zjnnNITZ1mRlZek73/mOxo0bp3feeUcnnniittpqK61Y\nsaLF63TnQH2JMWIAgEy148lSVp9mK2NSdj9p68nB4j43SWOmSLFcybKkvttIe10njfhmsH3UJMla\naJPwOmn4f3W4tLPPPluzZ8/W7NmzZWb62c9+prPPPlvLli2TFLQO1Q/I/8lPfqLbb79dzz33nNxd\nS5cu1bx587TNNttov/320wUXXKDq6mq99dZbuvXWWzVt2rSG65x44omaOXOm/v73v2vq1KkN66dP\nn64rrrhCc+cGNy6sXr1aDzzwQNJ6jzzySH3wwQe68847VVtbq5qaGpWXl2vevHmKxWL69re/rZKS\nEm3YsEFz587VjBkzWn3//uUQI33++ef605/+pMsvv1xXXnmlJGndunWKxWIaOnSoEomEbr/9dr3z\nzpc3SYwYMUKffPKJampqGtatXbtWgwYNUk5OjmbPnq277ror6fVnzJih0tJSrV27Vu6uxx9/XHPn\nztW+++6rkSNHauLEiTrnnHNUVVUld9fHH3+sF154Iem106q9fZnd+SXGiAEAQi3+TVg+2/3pg91n\n7eD+0A7uLxznvnbR5vvV1bjXrHVPJJqu37jM/YkJwbiyfxS4P7R9MD5s0UPtqu0rX/mK//Of/2yy\n7rTTTvNjjz02uMzGjX7hhRf6dttt5/n5+b7LLrv49ddf37DvQw895OPGjfN4PO477rijP/XUU+7u\n/sknn/iRRx7pgwcP9h122MFvvvnmJtfYsGGDx+Nx33333Ter6c477/Tdd9/d8/PzfcyYMf6Tn/yk\nYVssFvP58+c32f+DDz7wI444wocNG+ZDhw71Qw45xOfMmePu7suWLfMjjzzS8/PzfZ999vHf/OY3\nbY4Ri8fjPmDAAB8xYoQfccQRDe+p3kUXXeSDBw/2YcOG+XnnnedFRUV+6623urv7pk2bGt73sGHD\n3N39gQce8IKCAs/Ly/PJkyf7GWec4dOmTWuxhgcffND3339/Hzx4sOfn5/u4ceN85syZDdvXrFnj\np556qm+zzTY+cOBA32uvvfzee+9Neu3mkuUTdWCMmHkXz6XSlczMe3J9AIDuY2ZK+jeheoUUy5Fy\n8tp/4tr10uJ/SMtelvqOlsaeKMU3H7MF1Ev2uxiub1ffJkEMAJARWg1iQDfqyiDGGDEAAICIEMQA\nAAAiQhADAACICEEMAAAgIgQxAACAiDCzPgD0YtXV1aqqqlI8Hldubm7U5XRKQUFBt896DrSkI494\nSobpKwCgF6qrq1NpaanKy8sVi8WUSCRUWFio4uJiZWVlRV0e0Ct1ZPoKWsQAoBcqLS1VRUWFamtr\nG9bVP4R68uTJUZUFoBnGiAFAL1NdXa3y8vLNnpVX//zA6urqJEcC6G4EMQDoZaqqqhSLtfzPeywW\nU1VVVTdXBCAZghgA9DLxeFyJRKLFbYlEQvF4vJsrApAMQQwAepnc3FwVFhYqJyenyfqcnBwVFhZm\n/N2TQG/CYH0A6IWKi4slqcldk+PHj29YD6BnYPoKAOjFetM8YkBP15HpKwhiAAAAXaAjQaxbxoiZ\nWczM/m1mD4fLg8zsKTObZ2ZPmll+d9QBAADQk3TXYP2zJM1ttHy+pGfcfWdJz0q6oJvqAAAA6DHS\nHsTMbBtpaxACAAAgAElEQVRJxZJuabT6aEkzwtczJB2T7joAAAB6mu5oEfujpF9KajzYa4S7V0qS\nu38maXg31AEAANCjpDWImdkRkird/U1JrQ1eY0Q+AADY4qR7HrH9JR1lZsWS+kqKm9kdkj4zsxHu\nXmlmIyV9nuwEJSUlDa+LiopUVFSU3ooBAABSUFZWprKysk6do9umrzCzgySd5+5HmdkfJK1w96vM\n7NeSBrn7+S0cw/QVAAAgI/TY6StacKWkw8xsnqRDwmUAAIAtChO6AgAAdIFMahEDAADY4hHEAAAA\nIkIQAwAAiAhBDAAAICIEMQAAgIgQxAAAACJCEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACI\nCEEMAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhBDAAAICIEMQAAgIgQxAAAACJCEAMAAIgIQQwAACAi\nBDEAAICIEMQAAAAiQhADAACICEEMAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhBDAAAICIEMQAAgIgQ\nxAAAACJCEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACICEEMAAAgIgQxAACAiBDEAAAAIkIQ\nAwAAiAhBDAAAICIEMQAAgIgQxAAAACJCEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACICEEM\nAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhBDAAAICIEMQAAgIikNYiZWa6Z/cvM3jCzd83sinD9JWb2\niZn9O/z6VjrrAAAA6InM3dN7AbN+7r7ezLIkvSzpPEmHSqpy92vbONbTXR8AAEBXMDO5u7XnmLR3\nTbr7+vBlbni9L8LldhUKAADQ26Q9iJlZzMzekPSZpDJ3nxtuOt3M3jSzW8wsP911AAAA9DTd0SKW\ncPc9JW0j6UAzO0jSDZK2c/c9FAS0VrsoAQAAeqPs7rqQu68xs8ckFbr78402/UXSI8mOKykpaXhd\nVFSkoqKidJUIAACQsrKyMpWVlXXqHGkdrG9mQyXVuPtqM+sr6UlJl0p6190/C/c5R9Le7j61heMZ\nrA8AADJCRwbrp7tFbJSkGWZmCrpB73D3f5rZTDPbQ1JC0gJJ09NcBwAAQI+T9ukrOoMWMQAAkCl6\n5PQVAAAAaBlBDAAAICIEMQAAgIgQxAAAACJCEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACI\nCEEMAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhBDAAAICIEMQAAgIgQxAAAACJCEAMAAIgIQQwAACAi\nBDEAAICIEMQAAAAiQhADAACICEEMAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhBDAAAICIEMQAAgIgQ\nxAAAACJCEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACICEEMAAAgIgQxAOhFqqurtXz5clVX\nV0ddCoAUZEddAACg8+rq6lRaWqry8nLFYjElEgkVFhaquLhYWVlZUZcHIAmCGAD0AqWlpaqoqFBt\nbW3DuoqKCknS5MmToyoLQBvomgSADFddXa3y8nLV1NQ0WV9TU6Py8nK6KYEejCAGABmuqqpKsVjL\n/5zHYjFVVVV1c0UAUkUQA4AMF4/HlUgkWtyWSCQUj8e7uSIAqSKIAUCGy83NVWFhoXJycpqsz8nJ\nUWFhoXJzcyOqDEBbGKwPAL1AcXGxJDW5a3L8+PEN6wH0TObuUdeQlJl5T64PAHqa6upqVVVVKR6P\n0xIGdDMzk7tbu47pyUGHIAYAADJFR4IYY8QAAAAiQhADAACICEEMAAAgIgQxAACAiBDEAAAAIkIQ\nAwAAiAhBDAAAICIEMQAAgIgQxAAAACKS1iBmZrlm9i8ze8PM3jWzK8L1g8zsKTObZ2ZPmll+OusA\nAADoidL+iCMz6+fu680sS9LLks6TdJSkFe7+BzP7taRB7n5+C8fyiCMAAJAReuQjjtx9ffgyN7ze\nF5KOljQjXD9D0jHprgMAAKCnSXsQM7OYmb0h6TNJZe4+V9IId6+UJHf/TNLwdNcBAADQ02Sn+wLu\nnpC0p5nlSXrSzIokNe9vpP8RAABscdIexOq5+xozK5VUKKnSzEa4e6WZjZT0ebLjSkpKGl4XFRWp\nqKgo3aUCAAC0qaysTGVlZZ06R1oH65vZUEk17r7azPpKelLSpZImSlrp7lcxWB8AAPQGHRmsn+4W\nsVGSZpiZKRiPdoe7/zMcM3afmf1Y0kJJx6W5DgAAgB4n7dNXdAYtYgAAIFP0yOkrAAAA0DKCGAAA\nQEQIYgAAABEhiAEAAESEIAYAABARghgAAEBECGIAAAARIYgBAABEhCAGAAAQEYIYAABARAhiAAAA\nESGIAQAARIQgBgAAEBGCGAAAQEQIYgAAABEhiAEAAESEIAYAABARghgAAEBECGIAAAARIYgBAABE\nhCAGAAAQEYIYAABARAhiAAAAESGIAQAARIQgBgAAEBGCGAAAQEQIYgAAABEhiAEAAESEIAYAABAR\nghgAAEBECGIAAAARIYgBAABEhCAGAAAQkexUdjKzLEkjGu/v7ovSVRQAAMCWoM0gZmZnSLpEUqWk\nRLjaJY1LY10AAAC9nrl76zuYfSRpH3df0T0lNbm2t1UfAABAT2BmcndrzzGpjBFbLGl1x0oCAABA\nMqmMEftYUpmZPSapun6lu1+btqoAAAC2AKkEsUXhV5/wCwAAAF2gzTFiDTua9XP39Wmup/k1GSMG\nAAAyQlrGiJnZBDObK+n9cPnrZnZDB2sEAABAKJXB+tdJOlzSCkly9zmSDkxnUQAAAFuClGbWd/fF\nzVbVpaEWAACALUoqg/UXm9l+ktzMciSdJem99JYFAADQ+6XSInaKpJ9L2lrSEkl7SDotnUUBAABs\nCVJpEdvZ3b/XeIWZ7S/p5fSUBAAAsGVIpUXs+hTXAQAAoB2StoiZ2QRJ+0kaZmbnNtqUJykr3YUB\nAAD0dq11TfaRNCDcJ95o/RpJ301nUQAAAFuCpEHM3Z+X9LyZbXD3PzTeZmZTJH2Y7uIAAAB6s1TG\niJ3QwroLuroQAACALU1rY8QmSSqWtLWZ/W+jTXmSatNdGAAAQG/X2hixpZLKJR0lqaLR+ipJ56Sz\nKAAAgC2BuXvrO5hlu3uHWsDMbBtJMyWNkJSQdLO7X29ml0j6maTPw10vdPcnWjje26oPAACgJzAz\nubu165hkQcfM7nP348zsbUmb7eTu41IoaKSkke7+ppkNUNCydrSk4yVVufu1bRxPEAMAABmhI0Gs\nta7Js8LvR3a0IHf/TNJn4eu1ZvaegkclSVK7CgUAAOhtWu2aNLMsSc+4+zc7fSGzsZLKJO0m6TxJ\nP5S0WsE4tPPcfXULx9AiBgAAMkJHWsRanb7C3eskJcwsv5OFDZD0gKSz3H2tpBskbefueyhoMWu1\nixIAAKA3SuWh32slvW1mT0taV7/S3c9M5QJmlq0ghN3h7rPCY5c12uUvkh5JdnxJSUnD66KiIhUV\nFaVyWQAAgLQqKytTWVlZp86Ryl2TP2hpvbvPSOkCZjMlLXf3cxutGxmOH5OZnSNpb3ef2sKxdE0C\nAICM0NWD9et9IukVd9/QgYL2l/Q9BS1qbyi4+/JCSVPNbA8FU1oskDS9vecGAADIdKm0iM2QNEHS\nSkkvSnpB0kvu/kXai6NFDAAAZIgunUeshZOPlvRdSb+QNNrdU2lN6xSCGAAAyBRp6Zo0s5MkHSBp\nd0nLJf1JQcsYAAAAOiGVrsnlkuZLulHSc+6+oBvqqr82LWIAACAjpK1r0sx2lXSgpP+StKOkee4+\nrUNVtgNBDAAAZIoun9A1PGmepDGSCiSNlZSv4G5HAAAAdEIqXZNvSXop/HrB3T/pjsLCa9MiBgAA\nMkJa75qMAkEMAABkirR0TQIAACA9CGIAAAARIYgBAABEJJW7Jv9gZnlmlmNm/zSzZeEkrwAAAOiE\nVFrEJrr7GklHKnhA9w6SfpnOogAAALYEqQSx+scgHSHpfndfncZ6AAAAthipPLj7UTN7X9IGSaea\n2TBJG9NbFgAAQO+X6iOOBkta7e51ZtZPUp67f5b24phHDAAAZIiOzCOWtEXMzA5292fN7NuNL9DI\ng+0vEQAAAPVa65o8SNKzkia3sM1FEAMAAOgUHnEEAADQBdLyiCMzOyucR8zM7BYz+7eZTex4mQAA\nAJBSm77ix+E8YhMlDZE0TdKVaa0KAABgC5BKEKtvYiuWNNPd3220DgAAAB2UShCrMLOnFASxJ80s\nLimR3rIAAAB6vzYH65tZTNIekj5291VmNkTS1u7+VtqLY7A+AADIEGkZrK9gqopdJJ0ZLveXtFU7\nawMAAEAzqQSxGyRNkHRiuFwl6c9pqwgAAGALkcqzJvdx973M7A1JcvcvzKxPmusCAADo9VJpEasx\nsywFXZQKH/rNYH0AAIBOSiWI/a+kf0gabma/k/SSpCvSWhUAAMAWIKVHHJnZVyUdomD+sH+6+3vp\nLiy8LndNAgCAjNCRuyZTDWJZkkao0Zgyd1/U7grbiSAGAAAyRUeCWJuD9c3sDEmXSKqUVKegVcwl\njetIkQAAAAikMqHrRwrunFzRPSU1uTYtYgAAICOka0LXxZJWd6wkAAAAJJO0RczMzg1f7ippZ0mP\nSaqu3+7u16a9OFrEAABAhujqMWLx8Pui8KtP+CWFc4oBAACg41IZIzbF3e9va1060CIG9AzV1dWq\nqqpSPB5Xbm5u1OUAQI+UlukrzOzf7r5XW+vSgSAGRKuurk6lpaUqLy9XLBZTIpFQYWGhiouLlZWV\nFXV5ANCjdGnXpJlNklQsaWsz+99Gm/Ik1XasRACZpLS0VBUVFaqt/fI/+YqKCknS5MmToyoLAHqN\n1u6aXCqpXNJGSRWNvh6WdHj6SwMQperqapWXl6umpqbJ+pqaGpWXl6u6ujrJkQCAVCVtEXP3OZLm\nmNld7l6TbD8AvVNVVZVisZb/Xy0Wi6mqqorxYgDQSW3OI0YIA7ZM8XhciUSixW2JRELxeLzFbQCA\n1KUyoSuALVBubq4KCwuVk5PTZH1OTo4KCwtpDQOALpA0iJnZHeH3s7qvHAA9SXFxscaPH6/s7Gz1\n6dNH2dnZGj9+vIqLi6MuDQB6hdZm1p8r6VBJj0sqUvCw7wbuvjLtxTF9BdAjMI8YALStq2fWv1HS\nPyVtp+BuycYn9nA9gF6qefgigAFA10tlQtf/c/dTu6me5temRQzoZkziCgAdk5aZ9cMTf13SAeHi\nC+7+VgfqazeCGND9HnnkEVVUVDSZPywnJ0fjx49nElcAaEVHglibd02a2ZmS/iZpePj1NzM7o2Ml\nAujJmMQVALpXa2PE6v1U0j7uvk6SzOwqSa9Kuj6dhQHofkziCgDdK5V5xExSXaPlOjW7gxJA78Ak\nrgDQvVJpEbtd0r/M7B/h8jGSbk1fSQCiUj+Ja7IxYrSGAUDXSnWw/l6S/itcfNHd30hrVV9el8H6\nQDfjrkkgg9WukzZWSluNkLL7R13NFidtd01GhSAGRIdJXIEMkqiV3r5cWniPFMsOlgtOkHa/OFhG\nt0jLXZOdYWbbmNmzZvaumb0d3oEpMxtkZk+Z2Twze9LM8tNZB4D2y83N1dChQwlhQCZ4+3Jp0X1S\nojpoFUtUB8tvXx51ZWhDWlvEzGykpJHu/qaZDVAwQ//Rkn4kaYW7/8HMfi1pkLuf38LxtIgBANCa\n2nVS6V5B+GouaytpUgXdlN0kXfOI9TezWPh6JzM7ysxyUjm5u3/m7m+Gr9dKek/SNgrC2IxwtxkK\nbgAAAADttbEyefejZQXb0WOl0jX5gqStzGxrSU9Jmibpr+29kJmNlbSHpNckjXD3SikIawomigUA\nAO211YhgTFhLvC7Yjh4rpXnE3H29pG9LusHdp0jatT0XCbslH5B0Vtgy1ry/kf5HAAA6Irt/MDA/\nq2/T9Vl9pTHH0y3Zw6VyK4WZ2QRJ35P0k3Bdyvewm1m2ghB2h7vPCldXmtkId68Mx5F9nuz4kpKS\nhtdFRUUqKipK9dIAAGwZdr84+L7o3qA70uukMcd9uR5pUVZWprKysk6do83B+mZ2kKTzJL3s7leZ\n2XaSznb3M1O6gNlMScvd/dxG666StDI8H4P1AQDdo7fPs9Xb318Pl9Z5xMysX9hF2Z6C9lcwxuxt\nBd2PLulCSbMl3SdpW0kLJR3n7qtaOJ4gBgDoPObZQjdISxALuyVvlTTA3ceY2dclTXf30zpeaorF\nEcQAAF1hziXBvFp1G75cl9U36L77+qXR1YVeJV0Tul4n6XBJKyTJ3edIOrD95QEAEIHadUFLWOMQ\nJgXLi+4NtgMRSWlmfXdf3GxVXRpqAQCg6zHPFnqwVDrGF5vZfpI8nMj1LAUTswIA0PMxzxZ6sFRa\nxE6R9HNJW0taomBS1p+nsygAALoM82yhB0vrsyY7i8H6AIAuUX/XZJN5to7nrkl0qXTdNfm/Laxe\nLam80QStaUEQAwB0KebZQhql667JrRR0R34Yfo1T8ODun5jZde2uEgCAqGT3lwZsRwhDj5FKi9hr\nkvZ397pwOVvSi5L+S9Lb7r5L2oqjRQwAAGSIdLWIDZI0oNFyf0mDw2BW3Z6LAQAA4EupjFD8g6Q3\nzaxMkimYzPUKM+sv6Zk01gYAANCrpXTXpJmNkvSNcPF1d1+a1qq+vC5dkwAAICOk7aHfZjZI0o4K\nBu5Lktz9hXZX2E4EMQAAkCk6EsTa7Jo0s58qmE1/G0lvStpX0quSDu5IkQAAAAikMlj/LEl7S1ro\n7t+UtKekVWmtCgAAYAuQShDb6O4bJcnMct39fUk7p7csAACA3i+VuyY/MbOBkh6S9LSZfSFpYXrL\nAgAA6P3a9axJMztIUr6kJ9x9U9qq+vJ6DNYHAAAZoUvvmjSzrSSdImkHSW9LutXdaztdZTsQxAAA\nQKbo6pn1Z0gqVBDCJkm6phO1AQAAoJnWWsTedvfdw9fZkma7+17dWhwtYgAAIEN0dYtYTf2L7u6S\nBAAA2BK01iJWJ2ld/aKkvpLWh6/d3fPSXhwtYkDvVrtO2lgpbTVCyu4fdTUA0CldOrO+u2d1viQA\n0OaBK1ErvX25tPAeKZYdLBecIO1+cbAMAFsI/sUD0HnJWraSBS5PSIsfkBLVwZckLbov+P71S7u/\nfgCISLvmEetudE0CPVxbLVtzLgkCVt2GL4+JbSV5jeR1m58vaytpUgXdlAAyUlcP1geA1r19eRC0\nEtVBq1iiOlh++/JgeeE9TUOYJCU2thzCJMmygpY1ANhCEMQAdEyyoFW3QVp0r7R2QfvHe3ld0L0J\nAFsIghiAjtlYmTxoWXivTyLZzDdZQRdlk1V9pTHH0y0JYItCEAPQMVuNSB60vE4aMDYYL5bVt+m2\nrL7S2KlSwfHBmLDs/sH3MccFY8sAYAvCXZMAOia7fxC0mg/Gz+obhKrs/l8Gq0X3Bq1kXvdl4Ipl\nS7v+innEAGzRuGsSQMfV3zXZJGgdv/l8YEzcCmAL0JG7JgliADqPoAUABDEAAICoMI8YAABABiGI\nAW2pXSet/Tj4DgBAF+KuSSAZHkwNAEgz/poAyTR+fA8PpgYApAFdk0BL2np8D92UAIAuQBADWtLW\n43t4MDUAoAsQxICWtPX4Hh5MDQDoAgQxoCX1j+9p6TmJPJgaANBFGKwPJNPacxIBAOgCzKwPtIXH\n9wAAUsAjjpA+hBEAAFrVkSBG1yRax6SmAACkDX9J0TomNQUAIG24axLJMakpAABpRRDrTbr64dRM\nagoAQFrRNdkbpGscF5OaAgCQVrSI9QaNx3HVrgu+L7ovWN8ZTGoKAEBaEcQyXbrHce1+cTCJadZW\nQfDK2opJTQEA6CLMI5bp1n4slU1uOXBl95eKHpEGbNf56zCPGAAArerIPGK0iHVEVw+K74zuGseV\n3T8IdIQwdKHq6motX75c1dXVUZcCAJFI62B9M7tV0pGSKt19XLjuEkk/k/R5uNuF7v5EOuvoMj1x\nctP6cVyL7mvaPZnVN+hCJDihB6qrq1NpaanKy8sVi8WUSCRUWFio4uJiZWVlRV0eAHSbdLeI3S7p\n8BbWX+vue4VfmRHCpPQNiu8sxnEhw5SWlqqiokK1tbXatGmTamtrVVFRodLS0qhLA4BuldYg5u4v\nSfqihU3t6j/tEXry5Kax7GCW+0kVwZiwSRXBMo8gQg9UXV2t8vJy1dTUNFlfU1Oj8vJyuikBbFGi\nGiN2upm9aWa3mFl+RDW0TyZMbso4LmSAqqoqxWIt/9MTi8VUVVXVzRUBQHSiCGI3SNrO3feQ9Jmk\nayOoof2Y3BToEvF4XIlEosVtiURC8Xi8mysCgOh0e9+Vuy9rtPgXSY+0tn9JSUnD66KiIhUVFaWl\nrjYxKB7oErm5uSosLFRFRUWT7smcnByNHz9eubm5EVYHAKkrKytTWVlZp86R9nnEzGyspEfcffdw\neaS7fxa+PkfS3u4+NcmxPWsesfq7JhfdG3RHel0ww3yUd00CGYi7JgH0Rh2ZRyytQczM7pJUJGmI\npEpJl0j6pqQ9JCUkLZA03d1bHGDV44JYPSY3BbpEdXW1qqqqFI/HaQkDkPF6XBDrrB4bxAAAAJph\nZn0AAIAMQhADAACICEEMAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhBDAAAICIEMQAAgIgQxAAAACJC\nEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACICEEMAAAgIgQxAACAiBDEAAAAIkIQAwAAiAhB\nDAAAICIEMQAAgIgQxAAAACJCEAMAAIgIQQwAACAiBDEAAICIEMQAAAAiQhADAACICEEMAAAgIgQx\nAACAiBDEgB6kurpay5cvV3V1ddSlAAC6QXbUBQCQ6urqVFpaqvLycsViMSUSCRUWFqq4uFhZWVlR\nlwcASBOCGNADlJaWqqKiQrW1tQ3rKioqJEmTJ0+OqiwAQJrRNQlErLq6WuXl5aqpqWmyvqamRuXl\n5XRTAkAvRhADIlZVVaVYrOX/FGOxmKqqqrq5IgBAdyGIARGLx+NKJBItbkskEorH491cEQCguxDE\ngIjl5uaqsLBQOTk5Tdbn5OSosLBQubm5EVUGAEg3BusDPUBxcbEkNblrcvz48Q3rAQC9k7l71DUk\nZWbek+sDulp1dbWqqqoUj8dpCQOADGNmcndr1zE9OegQxAAAQKboSBBjjBgAAEBECGIAAAARIYgB\nAABEhCAGAAAQEYIYAABARAhiAAAAESGIAQAARIQgBgAAEBGCGAAAQEQIYgAAABEhiAEAAESEIAYA\nABARghgAAEBECGIAAAARIYgBAABEJK1BzMxuNbNKM3ur0bpBZvaUmc0zsyfNLD+dNQAAAPRU6W4R\nu13S4c3WnS/pGXffWdKzki5Icw0AAAA9UlqDmLu/JOmLZquPljQjfD1D0jHprAEAAKCnimKM2HB3\nr5Qkd/9M0vAIagAAAIhcTxis71EXAAAAEIXsCK5ZaWYj3L3SzEZK+ry1nUtKShpeFxUVqaioKL3V\nAQAApKCsrExlZWWdOoe5p7dByszGSnrE3XcPl6+StNLdrzKzX0sa5O7nJznW010fAABAVzAzubu1\n65h0Bh0zu0tSkaQhkiolXSLpIUn3S9pW0kJJx7n7qiTHE8QAAEBG6HFBrLMIYgAAIFN0JIj1hMH6\nAAAAWySCGAAAQEQIYgAAABEhiAEAAESEIAYAABARghgAAEBECGIAAAARIYgBAABEhCAGAAAQEYIY\nAABARAhiAAAAESGIAQCA/9/e3cdYWtZnHP9eszNMwzIqsYBtqVgqllAgyk4JBmqXUCiZZkWNpZC+\naVslKkFjaCXVBPqPFZvYUhobC9Zgi7TWdAXJGCAxWwqFsjO8ShWoCKmKC2tKHYIeZnd+/WOehdnp\n7rCzOzv3ObvfTzI559znefnNPPucc+193+c5asQgJkmS1IhBTJIkqRGDmCRJUiMGMUmSpEYMYpIk\nSY0YxCRJkhoxiEmSJDViEJMkSWrEICZJktSIQUySJKkRg5gkSVIjBjFJkqRGDGKSJEmNGMQkSZIa\nMYhJkiQ1YhCTJElqxCAmSZLUiEFMkiSpEYOYJElSIwYxSZKkRgxikiRJjRjEJEmSGjGISZIkNWIQ\nkyRJasQgJkmS1IhBTJIkqRGDmCRJUiMGMUmSpEYMYpIkSY0YxCRJkhoxiEmSJDViEJMkSWrEICZJ\nktSIQUySJKkRg5gkSVIjBjFJkqRGDGLSy+j1emzdupVer9e6FEnSAWa4dQHSvuj1eszMzDA2Nsbo\n6OiKbnv79u1MTk4yNTXF0NAQc3NzjI+PMzExwZo1a1Z0X5Kkg5NBTANpNULS5OQk09PTbNu27cW2\n6elpADZs2LAi+5AkHdwcmtRAWhiSXnjhBbZt28b09DSTk5Mrsv1er8fU1BSzs7M7tc/OzjI1NeUw\npSRpRRjENHBWIyTNzMwwNLTr02NoaIiZmZl93ockSc2GJpM8AfwvMAfMVtWprWrRYNmTkLSv88XG\nxsaYm5vb5XNzc3OMjY3t0/YlSYK2PWJzwPqqepMhTMuxGiFpdHSU8fFxRkZGdmofGRlhfHx8xT8Y\nIEk6OLUMYmm8fw2o1QpJExMTrFu3juHhYQ455BCGh4dZt24dExMTK7J9SZJSVW12nDwOPAtsB/62\nqq7ZxTLVqj71t9W8tMT+vESGJOnAkYSqyrLWaRjEfqqqnkpyBHAbcHFV3bFoGYOYlmRIkiT1i70J\nYs0m61fVU93tM0k2AqcCdyxe7oorrnjx/vr161m/fv0qVahBMDo6agCTJDWxadMmNm3atE/baNIj\nluRQYKiqnkuyFrgV+NOqunXRcvaISZKkgTBIPWJHARuTVFfD9YtDmCRJ0oGu2RyxPWGPmCRJGhR7\n0yPm5SMkSZIaMYhJkiQ1YhCTJElqxCAmSZLUiEFMkiSpEYOYJElSIwYxSZKkRgxikiRJjRjEJEmS\nGjzQ1CsAAAmQSURBVDGISZIkNWIQkyRJasQgJkmS1IhBTJIkqRGDmCRJUiMGMUmSpEYMYpIkSY0Y\nxCRJkhoxiGm/6PV6bN26lV6v17oUSZL61nDrAnRg2b59O5OTk0xNTTE0NMTc3Bzj4+NMTEywZs2a\n1uVJktRXDGJaUZOTk0xPT7Nt27YX26anpwHYsGFDq7IkSepLDk1qxfR6Paamppidnd2pfXZ2lqmp\nKYcpJUlaxCCmFTMzM8PQ0K7/SQ0NDTEzM7PKFUmS1N8MYloxY2NjzM3N7fK5ubk5xsbGVrkiSZL6\nm0FMK2Z0dJTx8XFGRkZ2ah8ZGWF8fJzR0dFGlUmS1J+crK8VNTExAbDTpybXrVv3YrskSXpJqqp1\nDbuVpPq5Pu1er9djZmaGsbExe8IkSQeFJFRVlrVOPwcdg5gkSRoUexPEnCMmSZLUiEFMkiSpEYOY\nJElSIwYxSZKkRgxikiRJjRjEJEmSGjGISZIkNWIQkyRJasQgJkmS1IhBTJIkqRGDmCRJUiMGMUmS\npEYMYpIkSY0YxCRJkhoxiEmSJDViEJMkSWrEICZJktSIQUySJKkRg5gkSVIjBjFJkqRGDGKSJEmN\nGMQkSZIaMYhJkiQ1YhCTJElqxCAmSZLUiEFMkiSpkWZBLMm5Sb6Z5NEkH2lVhyRJUitNgliSIeCv\ngV8DfhG4MMnxLWrR/rNp06bWJWgveewGm8dvcHnsDj6tesROBR6rqierahb4R+C8RrVoP/EFZXB5\n7Aabx29weewOPq2C2M8A/73g8Xe6NkmSpIOGk/UlSZIaSVWt/k6T04Arqurc7vFlQFXVlYuWW/3i\nJEmS9lJVZTnLtwpia4BHgLOAp4B7gAur6hurXowkSVIjwy12WlXbk1wM3Mr88OhnDWGSJOlg06RH\nTJIkSQMwWT/J5Um+k+Te7ufc1jVpaV6sd7AleSLJA0nuS3JP63q0tCSfTbIlyYML2g5PcmuSR5Lc\nkuSVLWvUru3m2PmeNwCSHJ3ka0keTvJQkku69mWfe33fI5bkcmCmqj7Vuha9vO5ivY8yP//ve8Bm\n4IKq+mbTwrTHkjwOrKuq/2ldi15ekjOA54DPV9XJXduVwA+q6pPdf4YOr6rLWtap/283x873vAGQ\n5DXAa6rq/iSHAdPMXw/13Szz3Ov7HrHOsj6BoKa8WO/gC4Pz2nDQq6o7gMWh+Tzguu7+dcDbVrUo\n7ZHdHDvwPa/vVdX3q+r+7v5zwDeAo9mLc29QXmwvTnJ/kmvtYu97Xqx38BVwW5LNSd7TuhjtlSOr\nagvMv2EARzauR8vje94ASfI64I3A3cBRyz33+iKIJbktyYMLfh7qbjcAnwaOrao3At8H7K6V9q/T\nq+oUYAL4QDd8osHW33NQtJDveQOkG5b8EvDBrmds8bn2sudek8tXLFZVZ+/hotcAX9mftWiffRd4\n7YLHR3dtGhBV9VR3+0ySjcwPN9/Rtiot05YkR1XVlm4uy9OtC9KeqapnFjz0Pa+PJRlmPoT9fVXd\n2DUv+9zrix6xpXS/yA7vAL7eqhbtkc3A65Mck+QQ4ALgpsY1aQ8lObT7Hx5J1gLn4Dk3CMLO84pu\nAt7V3f894MbFK6hv7HTsfM8bKH8H/GdVXbWgbdnn3iB8avLzzI+9zgFPABftGH9Vf+o+bn0VL12s\n9xONS9IeSvJzwEbmu9OHges9fv0tyReA9cCrgS3A5cCXgX8GfhZ4Eji/qp5tVaN2bTfH7kx8z+t7\nSU4HbgceYv71soA/Yf6bgr7IMs69vg9ikiRJB6q+H5qUJEk6UBnEJEmSGjGISZIkNWIQkyRJasQg\nJkmS1IhBTJIkqRGDmKR9kuSoJDckeaz7fsqbk7z+Zdb58+6rzK5M8pNJ7k4y3V2bR4skuTzJh1vX\nIWnl9cVXHEkaaBuBz1XVhQBJTgKOAv5riXXeAxxeVZXkAuDBqnrv/i9VkvqLPWKS9lqSM4EXquqa\nHW1V9VBV3dk9v6Pn64Ekv9G13QgcBkwn+WPgSuBtSe5NMprk7CT/nmQqyT8lObRb75Qkm7pet68m\nOaprvyTJw0nu765UvrjGoQV13J/kA137Wd0+H0hybZKRrv3bSa7oeugeSPKGrv0tSe7r1plOsjbJ\nryT5yoJ9XZ3kdxds5+PdOpu7+m/peg4vWrDOpUnu6Wq7fEH7R5M8kuR24BdW5ohJ6jf2iEnaFycC\n07t6Isk7gJOr6qQkRwKbk9xeVecl+WFVndIttwVYV1WXJHk18DHgrKr6URfUPpzkE8DVwFur6gdJ\nzgc+DvwB8BHgdVU1m+QVuyjlvcAxXS2V5FVJRoHPAWdW1beSXAe8D/irbp2nq2pdkvcBl3bbuBR4\nf1Xd1YXDH3fLLvX1JE9U1ZuSfKrb35uBQ5n//sDPJDkbOK6qTk0S4KYkZwDPA+cDJwOHAPcCU0vs\nR9KAMohJ2l/OAG4AqKqnk2wCfgm4mZ2/oHqh04ATgDu7YDIC3MV8j9CJwG1d+xDwvW6dB4AvJPky\n89+xuNivAn9T3fe5VdWzSU4GHq+qb3XLXAe8n5eC2Mbudhp4e3f/TuAvklwP/EtVfXe+lCXt6C17\nCFhbVc8Dzyf5cRcazwHOTnJv9zdZCxwHvALYWFU9oJfkppfbkaTBZBCTtC8eBt65h8suTC2760UK\ncGtV/dZOjcmJwNeraleT+X8deAvwVuCjSU6sqrll1rNYr7vdTvc6WVVXJrm529+dSc4BtrHzFI+f\n2M125hbc3/F4uKvhzxYO7QIk+eAe1C/pAOAcMUl7raq+BhyS5A93tCU5qRte+zfgN7s5WkcAvwz8\nx47FdrPJu4HTk/x8t61DkxwHPAIckeS0rn04yQndOq+tqn8FLmO+J+mwRdu8DbgoyZpu3cO77R2T\n5Nhumd8BNi31uyY5tqoerqpPApuB44EngROSjCR5FXDWUttYuLnu9hbg95Os7fbx093f6nbm582N\nJhkDNuzhdiUNGHvEJO2rtwNXJbkM+BHwBPChqrojyZuZHzqcA/6oqp7p1tllj1hVbU3yLuCGbh5X\nAR+rqseSvBO4OskrgTXAXyZ5FPiHbpgvwFVV9cNFm70WeAPwYJIXgGuq6tNJ3g18qQtom4HPLFUb\n8KHMfzhhO/M9gV/t5qV9kfk5X99mfi7Xi7/OEn+zHcOktyU5HrirG+acAX67qu7rtvsgsAW4Z4lt\nSRpg6aZNSJIkaZU5NClJktSIQUySJKkRg5gkSVIjBjFJkqRGDGKSJEmNGMQkSZIaMYhJkiQ1YhCT\nJElq5P8AbEhwesk/tu4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc635688390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_recovered = np.linalg.pinv(P).dot(data_transformed.T).T\n",
    "\n",
    "x_mean, y_mean = np.mean(coffee_v_thesis, axis=0)\n",
    "\n",
    "# Add the means\n",
    "data_recovered[:,0] += x_mean\n",
    "data_recovered[:,1] += y_mean\n",
    "\n",
    "# Visualising the data\n",
    "plt.figure(figsize=(10,10))\n",
    "original_data_fig = plt.scatter(x=coffee_v_thesis[:,0], y=coffee_v_thesis[:,1], color=\"#7a7a7a\", s=40)\n",
    "recovered_data_fig = plt.scatter(x=data_recovered[:,0], y=data_recovered[:,1], color=\"#ffaf1f\", s=40)\n",
    "plt.ylabel(\"Pages of thesis written\")\n",
    "plt.xlabel(\"Coffees consumed\")\n",
    "plt.title(\"Coffee vs Thesis (mean adjusted)\")\n",
    "plt.legend([original_data_fig, recovered_data_fig], [\"Original Data Set\", \"Recovered Data Set\"])\n",
    "plt.savefig(\"recovered_data.jpg\", dpi=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have not completely recreated the data set, as we only selected one of the eigenvalues we have lost some of the information describing the data.  You can see that most of the data lines along a diagonal line, it is this line that is represented by the principal component."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#End of notebook"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}