PCA_week1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 1: Mean/Covariance of a data set and effect of linear transformation\n",
    "\n",
    "In this week, we are going to investigate how the mean and (co)variance of a dataset changes\n",
    "when we apply affine transformation to the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "1. Get Farmiliar with basic programming using Python and Numpy/Scipy.\n",
    "2. Learn to appreciate implementing\n",
    "   functions to compute statistics of dataset in vectorized way.\n",
    "3. Understand the effects of affine transformations on a dataset.\n",
    "4. Understand the importance of testing in programming for machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are a few links for your reference. You may want to refer back to them throughout the whole course.\n",
    "\n",
    "- If you are less comfortable with programming in Python, have a look at this Coursera course https://www.coursera.org/learn/python.\n",
    "- To learn more about using Scipy/Numpy, have a look at the [Getting Started Guide](https://scipy.org/getting-started.html). You should also refer to the numpy [documentation](https://docs.scipy.org/doc/) for references of available functions.\n",
    "\n",
    "- If you want to learn more about creating plots in Python, checkout the tutorials found on matplotlib's website \n",
    "https://matplotlib.org/tutorials/index.html. Once you are more familiar with plotting, check out this excellent blog post http://pbpython.com/effective-matplotlib.html.\n",
    "\n",
    "- There are more advanced libraries for interactive data visualization. For example, [bqplot](https://github.com/bloomberg/bqplot) or [d3.js](https://d3js.org/). You may want to check out other libraries if you feel adventurous.\n",
    "\n",
    "- Although we use Jupyter notebook for these exercises, you may also want to check out [Jupyter Lab](https://github.com/jupyterlab/jupyterlab) when you want to work on your own projects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the packages that we will use for the week. Run the cell below to import the packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "from sklearn.datasets import fetch_lfw_people, fetch_mldata, fetch_olivetti_faces\n",
    "import time\n",
    "import timeit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we are going to retrieve Olivetti faces dataset.\n",
    "\n",
    "When working with some datasets, before digging into further analysis, it is almost always\n",
    "useful to do a few things to understand your dataset. First of all, answer the following\n",
    "set of questions:\n",
    "\n",
    "1. What is the size of your dataset?\n",
    "2. What is the dimensionality of your data?\n",
    "\n",
    "The dataset we have are usually stored as 2D matrices, then it would be really important\n",
    "to know which dimension represents the dimension of the dataset, and which represents\n",
    "the data points in the dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of the faces dataset: (400, 4096)\n",
      "400 data points\n"
     ]
    }
   ],
   "source": [
    "image_shape = (64, 64)\n",
    "# Load faces data\n",
    "dataset = fetch_olivetti_faces()\n",
    "faces = dataset.data\n",
    "\n",
    "print('Shape of the faces dataset: {}'.format(faces.shape))\n",
    "print('{} data points'.format(faces.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When your dataset are images, it's a really good idea to see what they look like.\n",
    "\n",
    "One very\n",
    "convenient tool in Jupyter is the `interact` widget, which we use to visualize the images (faces). For more information on how to use interact, have a look at the documentation [here](http://ipywidgets.readthedocs.io/en/stable/examples/Using%20Interact.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD3CAYAAAAwh5neAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnWuMXdWV5/+rXq7yowqXbbCNCwPBEMCiwUKJCUlEQ54k\nachklGdHnhHSfJlpJUpPBqdHGqmlRIPyoTWtTGs0rSSNJw3uPDoZII9u0SSdSaMOAQcnsTFQYGz8\nLttVLhvb9bq150Pdc/jvVXcv7ypX3YKT9ZMsn3vPOfvss8/Zdf9rr7XXlhACHMepJi0LXQHHceYP\n7+COU2G8gztOhfEO7jgVxju441QY7+COU2EuqoOLyAdE5AUReUlEts5VpRzHmRtktn5wEWkF8CKA\n9wI4COBpAJ8KITwHAMPDw+5gd5wm0tPTI/q7i/kFfxuAl0IIe0MIYwD+DsA9F1Ge4zhzzMV08MsB\nHKDPB+vfOY7zBqGtGRfp7+9vxmUc5/eODRs2mPsvpoMfAtBHn9fVv2vI17/+9ejz4sWLy+2lS5dG\n+9auXVtuX3bZZeV2T09PdFx7e3u5LRKbHzy20NraCgDo6OjAyMhIdFyxDwBaWmJBMzk52bC8tra4\n2fS1GS6zKE9EEEKIyud74TrpMjo7O6N9HR0dDY/TYyv8ubjuyZMnsWLFiujafNzExERUBtd3bGzs\nguU3qsf4+Dg0k5OT09o+F932XE6tViu3dZvyPl0PLlPXN3Wf+r3i80ZHR6N9Z8+eLbcPHHhdBA8M\nDAAAPvGJT+Db3/42jh07Vu577bXXojKK623fvn3afTAXI9GfBrBBRK4SkQ4AnwTw6EWU5zjOHDPr\nX/AQwoSI/CcA/wigFcA3Qwi756xmjuNcNBdlg4cQfgzgx1kXUrKWZeGyZcuifYsWLWq4rWUWSylL\n4uVKaH0cX28uptUW5U1OTqK1tTW6tmUqpOqkj+X6s3zUFObRyZMnsXjx4uhYS9ayZOfnAsQSldtq\nPqYj55pVfJxuj0bvTvFcuA30e8uw1J7JM2NWrFhRbh85cqTcrtVqZpnafErhkWyOU2G8gztOhWmK\nmwyYLpFY4mkZlJLlluy0pKAlkViqWVI+dY6+tt7Hddaj3HzffJ6+Fz5O7+Pzurq6ym09SswSmus0\nOTkZlWnJQkuupswDXQ+WltoLknqG+nur/FRbWffFZdRqteh90eVbzzqFPo7f7xMnTpTb7E1aunQp\nTp06VX7Wbc/mgYX/gjtOhfEO7jgVxju441SYptngVvQXuwoAoLu7u+F5lg1u2cXFNkd9NTrOKpO3\nLZtel9coYqogZZ/r8q3rafsxVQ+OgNO2aSrqim1kXUd9Xykb13oueiwgZdPmRqvp8nPR7kvrPWNS\nEYCAPaaSGn/i8nTb63eXIzwt/BfccSqMd3DHqTBNk+jnzp2LPi9fvrzcZvcOEMsTdqtomW9F86TO\n0/LLitxiLHeaJRktV0pKTmp5xuj6W9dmuA14osj4+HhD06VR/SwTgGV+rsQtzIGOjg6Mjo4m79uS\nv5qU2TMTU8G6ltXGuXAZ7BrjySW6HmfOnIk+64lXKfwX3HEqjHdwx6kwTZPoel4tj5zzqDkQyyJL\nhltyieVZKmLsQuWlIrf0XGhrNLWRTAwhTKsHf86V2vq8lOxsVK8crLbS+7gN+Fnr554a9a/Vaklp\nn+vp0Mdapo5FrrTnd3Mm7ZsyIzhHgp4EpM3HoaGhrGv5L7jjVBjv4I5TYbyDO06FWTAbXOcWY9i2\n4W1tm1puLbZVi5lUXV1d03JncVSRjsZqlNdN10nXy3Lp6LGAVA41y9Wm93G7WjOm+Dg9Q4/bhMvX\n5fF9a3uZI+Bmm+QhlcvNclVZ74BlI7N9XjyXWq2Gtra2qK30mEduYg4m19XG7mLtOtb14JlmFv4L\n7jgVxju441SYpkl0LWtZglgpeC3JmEoZDMSScXBwEMBU9Jw2FVJRXEB+dBZfS0vGRnJVRDAxMZFM\nm6zvpVGqYav84hop+Lrj4+PT2j9Vj9wortzkG1yPmbjkuHxd95SLyzKrCvkrIhgZGYneCW3SsVTm\na1lyPddk0eWx+ajroVMxJ8vMOspxnDcl3sEdp8J4B3ecCtM0G1wvvWIlU0zZetrOsWys06dPl9ts\nv+jjUskIdR0tV5g1461RGGtvby+Ghoam5RZPlW/Z/1wvvpZ1L1zf0dHRZHvr8YTUmIGuh/V9bq5v\nyyVnuQa5zlbShUZhrN3d3Thz5ozpEuXny7a6fpY8NnD+/PloX+p56vvnMvVYA7/fFv4L7jgV5oId\nXES+KSIDIrKLvusVkcdFpL/+/3KrDMdxFoYcif4ggP8J4P/Qd1sBPBFCeEBEttY/3z+TC+eu4Ggl\nWmDZYq2maUkdPk7Xg5NUaFmbQkv0JUuWlNuFpOvt7cXp06cjiWe1B8tJXf9U++i24jJ13vnU7CaN\nlS88FVGmnwubRMPDwwCmVpMdHByMorOsvN/W0kJcD2437Q7l9imSLnR3d+Ps2bOmW4v3cTSmrge3\nqZXTPCX5Ozo6ometZ49pt1mKC/6ChxD+H4BB9fU9ALbVt7cBuDfrao7jNJXZ2uCXhRCKldKOAshL\n8eg4TlORnCgbEbkSwA9DCBvrn0+FEC6h/UMhhMgOHx4eLgvu7++fq/o6jkNs2LCh3O7p6ZkWDjhb\nN9kxEVkTQjgiImsADFzohK985SvR54997GPlthUuyjaKNctKuzMarW+2du1avPTSS8lraffFyZMn\ny222u7WNzPXQs+S4HsUyyW9961vx/PPPR/biJZdc0vAcwF4TjI/la2sbnMsotk+fPo3u7u7k+ILl\ngtJwm6SWEgZet7uB19fluv7667Fnz57IzmQ3kE7Yyej2TmXv4bEQoPHy1cVz4Ta1Ek9aa8ZxGfr9\n5gSKjdyv1157LV588UXs2lWOa+OJJ56Iyjhw4AAAYM+ePbCYrUR/FMCW+vYWAI/MshzHceaRHDfZ\ndgD/CuA6ETkoIvcBeADAe0WkH8B76p8dx3mDcUGJHkL4VGLXXTO5UK4MB2KpzPLRmjlkuXesifRc\npq4H5263kj/muix0pFJuwgcrkWBqKV3rOL2sU2oml75nK6KO93Gb6mdmJV3k9uB84frd4cgw/TxZ\nsqfy6wOxfGc5vWjRItMk4jpa98loky41i45NoLa2tujedBtkz1DLOspxnDcl3sEdp8I0bbKJNXFB\nS7/URAMtS3hkVEt03sdyb926dclr6Qi11Mi5HrlNjSADsQxlOdbT0zNNGqbuhaWmlqRcpo6EStVD\n5y1PmTpWogkr5zhfW0tjPo+X3+np6UmugKrhe9H32SjXWqN6MPyu9PT0RG2s5Tq3CY/068lU/O5Y\nE1as1UX5s15N9ODBg4m7ifFfcMepMN7BHafCeAd3nArTNBuc7WDAjgJKubx0GWyjXHrppdG+IjqJ\nrzU4OBit/wQgmRNc15Gjj/RMJ54FpfNVc5mFjdnX14djx46ht7cXjbCSRmjblz/zvek2TCW9yJ05\nBaTHE/Rntov12Au7uNiebW9vj+xWPk6PjeS6FNmW1svt8phHynUHTB9vSbWxbm+O2LMSN1rtz++H\njsRLvTsa/wV3nArjHdxxKkzTJDovFwzEkkbLOJadLP20vGYpqJcgZgnGkxi0e4Gl1PHjx6N9R48e\nLbdfeeWVcpsnhgDAsWPHym2W8gCwcuXKafW99dZb8cwzz0Rtwu67VatWRWWwpLYkuiU1rUi2VD61\n3FzwQPwMrXx7fF4x2WT16tU4ceJE1P7sgtImEbfPvn37on36GRZo04nNveLd+fznP4+HHnoIN954\nY7lv9erV0Xlr166ddh4wXcpz/bVrk0ktS1Wr1aJ3R0t069kw/gvuOBXGO7jjVJimSXSeuAHEElrL\nTi3FC7RM4eP0SOWrr75abhdzZjdv3oynnnoqOo5H31mGA7H841Hda665JjqOZb6O/tq4cWO5zSbF\nxo0bsXfv3obX1pL02muvRYrU3Gst4bQs5+3c1Tr5PCv9tFUGt1VhLm3cuBEHDx6M5Ds/2/Xr10dl\ncJIDDV/7iiuuKLe3b98eHcdmFick6e/vjyT6jh07ovN+97vfldu33npruc1yWtdfv5spL5E2c1L5\nAgA7Mi+6VtZRjuO8KfEO7jgVxju441SYptngHFkGxPaGnhGUWhLGsjv2798ffWZbqXCdbN68GZdf\nfnl03E033VRuv+Utb4n2PfLI65mo2GX25JNPRsexK0WPH/DYA1/72muvjZIEHDp0qNzW98l2q7b1\n2F5nd4yV162wdVtaWjA6OppcElc/F54xpW18PZsqBee507PpePyCZ0/p8RtuY+1SZHuaxzj0c/ng\nBz9Ybu/evbvcvu2223D77beXn19++eXovF/84hfl9q9+9atye9OmTdFxHDlnLXudimRraWmJztM2\nuL7vFP4L7jgVxju441SYpkl0HezPslDLyZRE18dxhJqOVOLIIpbQWhJdd9115baezHLVVVeV21/9\n6lfLbS33uF7WpJfBwcHymMHBwUh6pvLQAa9HfAHTIwJZ5nNklZWjrnDbLF68GCMjI1FbsSvSkvkD\nA3Gm7NRkFu025Agvds/pHHVslujJNxwteMMNNyTrwa7Svr6+6Dh2oXEU4fvf/35zUg2/E1wPNrGA\n+LnrdkwtpaVztfHnNWvWRGXoyLkU/gvuOBXGO7jjVBjv4I5TYZpmg6dmLAHT7cVUzmo9K4dtU21b\nszuJr10s+VLw2GOPldt33HFHtI/tQA5P1Tmw2R7SriW2CXkJntOnT0chuuwGsXJss5sJiO1Trode\n7ie1BFStVouul8rVruuoZ5PxGIiVHIPHCfjZrl69OmoPLkOPSfCMQJ34gJ872//63WH3V+FyWrVq\nFV577TUcOXKk3KevnVpiykrYqcd2UmsE8LXa2tqi+ms3mU4emsJ/wR2nwuQsXdQnIj8TkedEZLeI\nfK7+fa+IPC4i/fX/l1+oLMdxmkuORJ8A8KchhF+LyDIAO0TkcQD/DsATIYQHRGQrgK0A7k9eSLk6\nWIanls7R+3hGFxDLIO2KYEnG0lXPzGIZ/pOf/CRZD3YfabnEslObIizRuR6dnZ1Jea2jwljm6jZg\nd4wlr62ED6lkDTp5RSqHt64/o6Ur11+7L7mNWZ5qmcwSt3A9Nroel6ffMTZL+JzR0dHIDaoj9vi+\n+Z3TpoiVU57rkloNtb293cydd/311yfLZy74Cx5COBJC+HV9+wyAPQAuB3APgG31w7YBuDfrio7j\nNI0Z2eAiciWAWwA8BeCyEEIxGnEUwGWJ0xzHWSAkd5VCEVkK4OcAvhJC+L6InAohXEL7h0IIpR0+\nPDxcFswTABzHmTs4+UVPT880WzfLTSYi7QD+HsBDIYTv178+JiJrQghHRGQNgIF0CfHMHiAOudRu\nhFS+aW1/ctijziLCNmFhU4nItDBTtp20zZlaX0rbZWxv6TpyvQoba8OGDejv709mSNEZQPh6eiyD\n245DX7UrhtuxGCfo7u7G6dOno7biGW88PqGvrevIyQ7ZlcfPCIjDi4u2ufPOO/HTn/40GtvgOulx\nDf5srePGNrKVa764bldXF86fPx89M+2W5GdtrZun16FjuF6NsvAsW7YMZ86cid5N/c4dPnwYgJ3d\nBsgbRRcA3wCwJ4TwF7TrUQBb6ttbADyiz3UcZ2HJ+QW/HcBnAfxORHbWv/szAA8A+I6I3AdgP4CP\nz08VHceZLRfs4CGEfwGQ8mPdlX0hw02mYTnCUlBHeFluG5aohbvk3Llz01xcbCpoOcaSnSPDtAzn\n+uqIKZZ0emYcSzC+T91Wqdl1QCxluQ20dNUzlXibJa+W3gzv08+vUXsDdgQjuwNbWlqS0lvfS8q1\npOuRWk4JSM/2WrRoUXSeNh+tKEPGunbq3df3zKaCvq5+D1J4JJvjVBjv4I5TYZo22USPclsrRFp5\ntVPH6VFM/swyWUeJpY4DYinFx2mJlYoE0+exzBobG4uOTUXeAfmjwVaeO743a1kjHjnX92LlT0+Z\nS/o4lrx8XytWrIi8FtZEEW5TXX4q0YV+to0SLSxbtgxnz541k2WkTAc9yp2bXz4V0akluo7m0+9B\nsr5ZRzmO86bEO7jjVBjv4I5TYZpmg88EtkXYDrFcPxrex3nAdXSWNesnNZvHsq+s5BV8LW0jW3ar\nZYOn7H9tZ1shydYsNIbdZJabhvfp9eTYzuZ76erqisYeUlGEQNx2ur1Tz0a3B7s9i2dURJDxs0gl\nZ9DXtsZlLPs8NRY1MTFhupJzQ8z9F9xxKox3cMepMAsm0a0IJ2YmcjhFIfe6urqmyT0rB3YqskrL\nZC2tGK6/lqS8j2Whll98nuUeSS0RrLEipqx89Yy1jJQVYcjSW5sXLIe5DOte9PvB17PceimzRERM\nF26qXlb5ltlm3Vsq+rDR5xT+C+44FcY7uONUGO/gjlNh3hB50S37gkMKtT1nuclS9qh2i1k2LV/P\nWjvMsn1TLi6deCJ1LcB2x+SOQ+gwyEbX1eh7sVxGqbEMbcfrfO18Tuo8fc/cPtY7wOfp8F9+hjoJ\nouU2TLlwNdY+LoPbTb9H1juXG87tv+COU2G8gztOhVkwN5nlKkjJLi33rLzRqcnyugzep6Ugy/nc\nPNeWuaFdJ6nzLCmsTYpUvWYr5Vky5spfXS9+LtoEYNOEZ/ZNTk5GZeRKXC1d+XOjPHSNysh9LkB6\nJpiurxXJlhs5mOtKtvBfcMepMN7BHafCLNgouiUZU5FbWp7mTsxPTfgAYjmsl9nhMjjpgCXXrVFp\nvpY12USXwTJO1zEl47QsbJRLrKOjAyMjI9H1rOV+rNH8lNS0pDynE+7p6cmWpJaU59Fyfk7WZBC+\nL/1crMQWVqQZY72nKZnf2tqaHSFp4b/gjlNhvIM7ToXxDu44FaZpNniuKwyI7SA+L9eW0efxtuWC\nsmZPcX11GbzPSkpR2G8hhGkRU2w76nqcOnWq3NbL5aZcOtpGY5cU29yjo6PJ2V+5s+usa2vbkZNG\n9Pb2lttLlixpuIwPMH38wEqAmZrtZR2XG9mny8ltq9wEjJpUlJveZ+G/4I5TYXLWJusUkV+JyG9E\nZLeI/Hn9+14ReVxE+uv/L79QWY7jNJcciT4K4M4Qwmv1VUb/RUR+AuDfAHgihPCAiGwFsBXA/alC\nLOmqJQxLYD7OikayyI1MsuRZ7mQZK5KtkIXj4+Noa2uL7pMlupa1J06cKLf1CqgpeaknV/BnnT/9\nQrK0wHoWqQQQuaZTCCFy0Wl3YM619PW4bazEEIyV1x5Iu8k0VrRao3yBmsnJyeg4XV9riamoHhc6\nIExRGHDt9X8BwD0AttW/3wbg3qwrOo7TNCQnO6OItALYAeAaAH8VQrhfRE6FEC6p7xcAQ8VnABge\nHi4L7u/vn/OKO44Trw/e09MzTdJmjaKHEGoAbhaRSwD8QEQ2qv1BRMy/FPv27Ys+8+qOesVP/mzJ\nwtz54Jw22VriyJobbkl0a7SW9xWj4+Pj42hvb8+W6K+88kq5ffr06WhfrkRvlD+sp6cHw8PD0QQQ\nji6zJt/o9khF9+njBgYGyu1iFL2rqwvnz59PSnQ9ip6ayw3EbcDH5Uj09vb2afeh35fUvP3cyD59\nXiOpXaRv5tVudb2K87iDN2JGbrIQwikR+RmADwA4JiJrQghHRGQNgAHr3NzZTED+pHpGN0CjhIa1\nWi3bvZNbP31ta1ZRsT0+Po7W1tbkjCPdifklsHKEp1w/us7aLcbrhaVsWL3PeoGtDtjd3d3wOJ0H\n3Mp9bi1jnBpTyV3DrKWlxRwfSqGfy2zscy5jfHzcbGMrXDq61oUOEJFV9V9uiEgXgPcCeB7AowC2\n1A/bAuCRrCs6jtM0cn7B1wDYVrfDWwB8J4TwQxH5VwDfEZH7AOwH8PF5rKfjOLPggh08hPBbALc0\n+P4kgLtyL6Qli5VYgKWKZZvmunf4mlYkmxUtZNlz1ow3K7EAl892pnadcL10LjeOerMiq1Kzp7Qk\nZbTtmzt7KpVzTF9bf8/X49l7udfSdbaWdc4tz3LvWjKcn4WW0/wu8fty9uzZcnt0dDQp3/V5Fh7J\n5jgVxju441SYpk02mcmILMs6a4Sa0dK70Qjk2NjYtDKs/F6p4yxy3SOTk5NJCWaZItos4RFwKyED\nf9YmQCoRhZbo1gh7asVM7QZKTcyZnJyMrpeacKTP08+Ty7A8N3xecV9FXriUxwFIj15b75XelzsC\nzm1ltYGF/4I7ToXxDu44FcY7uONUmKbZ4JYdYtm+1sR5to+0OyYVsjgT11oqjNVaukiTmrRvhd1q\n+5ntbHYfAdNDUlOwO61wx9RqNSxZsiQKEeVrWQkedf3Z9rWSdHDb6bEWPpbfD/3MUkkXgNiNmEog\nodGzAXOXTUqVobHCnFNRbSEEMwQ6F/8Fd5wK4x3ccSrMgrnJrJUTZxPsbyWNsFa7tBILpOSTlRRA\ny7FUgoNarRbVy5LoPEFj+fI4cc758+fLbY6E0vXgWWKFzD927BhWrlwZJZFgaazlKLedXiWUy+Dj\ntAnBbcplhBAilxo/P/3MrHcnlXzDivwqjqvVamhtbc06FrCXJ7LMFO0ebHSOluj6PnNNM/8Fd5wK\n4x3ccSqMd3DHqTBNs8Gt0ElrCVi2NSxXgQ5VZTsnlSv7QvtS7q/c7C5AOjTTCj3UCR9WrFhRbusx\ng9/+9rflNtvB69evj47j9lm9enX0Pdu41iw/zjDC20CcGJLzuL/rXe+KjuNnODw8XG6Pj48n22q2\nSyFbYzmNwl3HxsawaNGiqA20PZ5a+yw3CaWuVypjUa1WM13EnBHJwn/BHafCeAd3nArTNIn+85//\nPPp82223ldtWtBDLFB21ZMnrlAlguTNyky7q46wcZEwhC4vccI1mNAGxdAWAJ598stxmea2vvXLl\nymQ9WPbzdY8fPx5JdGsmGJfBSyHpenBixf3790fHcdRcsQzT6tWrMTg4GLkA+VlbM7qs9rZcrNY+\nK8+9tUxVCivpBd8LRwOOjY2ZkXg6P34K/wV3nArjHdxxKkzTJPrDDz8cfX7HO96RPJYlDEs/HdFk\nSanUCKSVeMKCZaKesG9FTDWSgoVET5kE69atiz5zTnkdQcaTQ3hbl82SuqhvZ2cnTp06FeWhZ1nI\nkXFAWuYDsfS+8sory+2jR49Gx/FkGX6ebW1tUT2s/HK5S0VZ0WQpb0GRsz5FzkIhunxrH3tFdCQb\n38szzzwTlbF7924AwN13323Ww3/BHafCeAd3nArjHdxxKkzTbPDDhw9Hnw8dOlRuX3fdddG+2SZo\nYFJL2Mw24YMV4cVYUVd6pha7odim1ckOeTaZzovOEWp6WWAmtaxRe3t79JltdW2Dp5IzArFLzXI9\nsg3OUWGrVq2K7jM1G1CXb+2zItka2eojIyNob2833aVWkhHGcqelxgl0Hnd2lz733HNRGXv37k1e\nm/FfcMepMNkdXERaReRZEflh/XOviDwuIv31/5dfqAzHcZrLTCT65wDsAVDoqK0AngghPCAiW+uf\n70+dzAkHAODHP/5xuX3jjTdG+1j6zGRVUiYl8XITSOh6sNlg1cnK16Yj2dhFYuU7Y6wEFamINCC+\n70Iaj42NTcvxxuex6wuwl+1JTe7Rrs1UIodFixYl3Y36Wmxm6WfB7ZNrSunVVtmUys3/rl2n1sQi\nrpc2xxiOCDxy5Ei0z5rMEtUj5yARWQfgQwC+Tl/fA2BbfXsbgHuzrug4TtOQHMe9iHwPwH8HsAzA\nfw4hfFhEToUQimWFBcBQ8RkAhoeHy4L7+/vnvOKO4wAbNmwot3t6eqZFbV1QoovIhwEMhBB2iMgd\njY4JIQQRMf9SvOc974k+cyTbF7/4xWgfSzceWdVy0op2aiTRX331VfT19VnVjNDB/wVWDrmcecBF\ntBTLSa7/0NBQdB7nXdMmhpa5jbZ1+SzROzo6IonKkxh01Jwl0flYvi9dDx45L57n4sWLce7cuSiS\njc0DS2pbEt2Kgmw0H3xwcBC9vb2mRE/lSdOmE+/TP6Lcxvysi/u8+uqrsXfvXjz//PPlvu3bt0dl\nFFGFenRdk2OD3w7gj0TkbgCdALpF5G8BHBORNSGEIyKyBsCAVYhOFrhr165yW3cKfoktezQ3ZNFa\nU8t6EIxecjenTkA6iaG+59msV6XP4z9qlk1YvIwiMu3FtNw77ArTzyX1nHQZqbERK2TYStKhyV1W\n13Kh8b3MNt85l2ktH5z6fmJiInKFaZflnK1NFkL4UghhXQjhSgCfBPDTEMIfA3gUwJb6YVsAPJJ1\nRcdxmsbF+MEfAPBeEekH8J76Z8dx3kDMKJIthPDPAP65vn0SwF3ZF1LSiXN6HTx4MNp39dVXl9tW\nsgZGy6yUO8xye+QuMWPN1LLkWLHd0dGB0dHRpPyzcrxbObYZa2mkwh4NIaClpcWcdZUq00qwoWeJ\nMak88SGE5BLKM3GVppJo6DZlc0NHN2p7neGxGGt5X762doWllnni9+i1117Dyy+/XH7W75WOJEzh\nkWyOU2G8gztOhWnaZBM9csguL3YHAGmJrmWKFbnFWMsCsTyzZJYVaWZNfphNQgktEfmzlVPOkoyN\nogM7Oztx/vz5SO5Zo+jcVlYEFqf01c+MP+tIMB7RbxR51+g8i9SkDiA9aaSIMizQ7Z07wp4yRRqV\nWVCYqps2bcLBgwejSDbLXWfhv+COU2G8gztOhfEO7jgVpmk2uHaXsN334osvRvve9773NSzDsj81\nqaSLuTaxxlpqyXLHNAqJnJycRGdnZ7ZLh9G2bypSTkc+cZ2LpYWuuuoqHDlyJJnIwXLF6HECfr5W\nsgauL5/T1tYW1dGy8S2XJT93q/5W21v5yFM2uJU/XUcLchncHi+99FK0zUku9X3mvsf+C+44FcY7\nuONUmKZJdCtRgY5kY/cA5yDLdTfoY3Mjjqz8XrkTHKzcX4VkPH78OFasWJGMirIityzJaCUZ4Gtx\nFNfZs2cjOc/SWyd8aDQjLfW5UXlAWsp3dnYml4fKzTEOpNtgJktWsalgRQRaCR9SphOQXqaKl3na\nv3+/mfhG7tbaAAAQNElEQVTEahPGf8Edp8J4B3ecCuMd3HEqTNNscG2bsl1SLCNbwGtx8TpXmtQa\nZjPBCilkuP5WthTLftb52VP1t5JS5CaN5HXKLLTtbGWFsZIkpMYo9PfcBuw21PnIrTGP3LBktoP1\nfaZs9RBCdgYZy1bnMQ89tsPnHThwoNzmsaeBgYGojfX7bc14i+qbdZTjOG9KvIM7ToVpmkS38miz\nqwCIZ5dt3ry53NaRSZzojxP2AelINo21JFHqPJ0gwIrAYpnI5509eza6H3ZdcQQTELsYOQGjPs9C\nn5eCy7PyxOnr8r1xDnwtTxvJ8NHRUbS3t0dSmdsmNy8fkHZdaTctz2bUiSZSS1Zp+H3Rx1mRj/yZ\nkzrwc9fvwGzxX3DHqTDewR2nwjRNolujy3pEkCefHD16tNzWkVUs97ScTOVM1yOwjSaDFKRGr/W1\nOBKsmMhRcPz4cWj6+vqwa9euaGSXzRRdhiU1UznTOZ+8LqNog2uuuQaHDx+etnRPwUwmOHBb8Wiw\nHs1nU6qQ8q2trTh37lzUHmwCWCu2WpGD1vuRyvE+OTlpvi+piSIabg9tznDuNW4rnXjCmrQzp0sX\nOY7z5sQ7uONUGO/gjlNhmmaDazcC207a/fXCCy+U27xwoY5qs6LQUm4yK2GitjlTNpx2YfA4gZ4Z\nx/sKW6yvrw87d+5Mumq0nZ1abwuIXYVc3xMnTkTHNbId3/3ud0dJBoB4nGMms7gYbm/tUmQ7e9Wq\nVQCmXKFPP/00VqxYUe6zcquzG06PNaTOs2Yi6vcolWxTf7aWvbKSb/B4C9vnenmsuUjm6b/gjlNh\nsn7BRWQfgDMAagAmQgi3ikgvgG8DuBLAPgAfDyEMpcpwHKf5zESi/2EIgXXfVgBPhBAeEJGt9c/3\np07OzeMMxLLuscceK7ff/va3R8exVLOC73OljrXPco+wHOvt7Y32sQxn90hfX18kz1juadcS19+S\n7yyNrYg6PqezszO5TLK+TytqjuvPMly3Kbdj4Q7cvHkz+vv7cejQoYZl8HMGgNWrVze8LgCsXLmy\n3NbmAZNKBjE5OZmd09xKyMBLBPNzB9LJN/S7aMlwa0IMczES/R4A2+rb2wDcexFlOY4zD0jOL6uI\nvAJgGFMS/X+HEP5aRE6FEC6p7xcAQ8VnABgeHi4L5oEyx3Hmjg0bNpTbPT09037ycyX6O0MIh0Tk\nUgCPi0i01lAIIYiI+Zfiox/9aHxhY74vjwxfccUV5faXv/zl6Li+vr5yW0e5NRrJPXDgANavX5+8\nrhUVxRFjPDIOxNFqWsbyeYVUu/nmm7Fz586kRNfMl0T/9Kc/jYcffjg61orYm2uJXjyzz372s/jW\nt74V5d/LlejFSHxBSqLr9uX3r5C7IyMj6OzszM61xu+pluj8TuhoxqGh14eqvvvd75bbhw8fBgD8\n6Ec/woc+9CEz3XdRZ73slyarg4cQDtX/HxCRHwB4G4BjIrImhHBERNYAGLDKmMkSsOw2e+6558rt\nX/7yl9Fxlr3FL4iVAM8Ke0y5grRbj9fistw2/JLddNNNUUgqb1s5wa0XjtEvNP8h4O1LL700ahPu\nZLqt2K7UbZPKd67h87htent7kx1cj0lYdUwlWshdFrlWq0VlWu3IfwiscFTL1ZaL9QfK4oK9TkSW\niMiyYhvA+wDsAvAogC31w7YAeCS/uo7jNIOcPwOXAfhB/a9hG4CHQwj/ICJPA/iOiNwHYD+Aj89f\nNR3HmQ0X7OAhhL0A/qDB9ycB3DXbC+cO87P99eCDD0b77rjjjnKb3VFALKNZnlkT83NdaFqG87W0\nvG40O2toaAgrVqyIIrdSch2I7XhLkrKMs9qXj1u7dm1kS1rJMbiNtWnAZXJ7WFKSZfj69eujcRTe\np12gvC83d9tMcurnJvBgdNISfmZWHa2IulSe+JngkWyOU2G8gztOhfEO7jgV5g2R0UXbxWxvsO2r\nZ0i9+uqr5Ta7ToB0ZhJrbS8r2wvbgdpOZZtQzxyysn6wrcr+XO0W4jJ1W7H9zMfp66ZcM0uXLo3a\nKtc3r8c8uE34WWiXIl+Lj+vr64vaOHfpXysc1fKDc3sUdraIYGxsLHKFaRs8lQvdSpKo3xf+zO+j\nDhnm9tH9pxmhqo7jvMHxDu44FeYNIdE1LH0sCcaRVVq+s8RjOWlFsuVGuVnJK7SpwBJal59aMljf\nM7uPtCuP740j6nRkFdeZ9y1btiwqP7WkMZDO8a7rlTKxgLT7q6urK9pnJeJIXUvX0TK/mEJ2t7e3\nY2JiInJxWe5ANom0lOf3zwqBTi0H1draauZWn7NINsdx3rx4B3ecCtM0ia4lRSovNRDLKZauenTZ\nilDjFUs5YkzPxrJGWlOy3Mr/Zq2myfepV9O0Vu7MHem3ZnGlcrwvWbIkOWtpJvVgKcvb2jTjOmqz\nJDW6bC0VpSU0S2UrJxufV5gs7e3tGBkZMctn+c4Rh/q5873pfaloQS3RGcvTZOG/4I5TYbyDO06F\n8Q7uOBWmaTa4xrJRUraptkN4Ur22sdhWKuxuEZlmg1vL1LIbx1ofjOurZz6l7K2WlpboPCsJgFVH\ny63F6CVyi3InJiaiMq31sFLJK/SxfM/aBm8UoVar1dDW1hYda83y42et3VOpbCxWBGBRxrJlyzAy\nMhKVqZdd5neOy9cZhZjc91vfM7exHofwSDbHcbyDO06VWbBINpYfViC95S6x8qnxeUXEW3d3t+km\ns8pP1V0fp02FVNKBEEIyQYMuw5K8LJs5Qs1yp1n5t61c34yVHCOVbONCpGSnluF8n1p65y6ry9Kb\n20pHsmmJztGTnIyE3X9AbMZZ8jrl6tX3pdvRly5yHMc7uONUGe/gjlNhmmaDWzOpcrHsOcvOYXtI\nT8y33HUp95E1u8lC275cvtUes8mjbWEtl2slbrRCKVPhxXosoJE7rVarobW1NZm/PqeMglQZOhFH\naq250dHRqA20Dc7X45l8+t20lsfmY1PjJiGE6DwdMmvNsoyulXWU4zhvSryDO06FWTCJbkljlki5\nOaW1DGJZbkkdnhGk852nlr6xXH6W+0i7/1LRWlb+MMvNZ2GZANwm1tpkXA99nyl3T67ZMz4+nnRx\n5Ubo6c9WbrVUm46MjETRarq92TVmzVazXIV87dS7o9vXShph4b/gjlNhsjq4iFwiIt8TkedFZI+I\n3CYivSLyuIj01/9fPt+VdRxnZuRK9L8E8A8hhH8rIh0AFgP4MwBPhBAeEJGtALYCuD9VgBWtZu2z\nJmRYub/0pI8CnTTi5MmTDa8LxJKdJZiOWrJGjXmUlGWWFXGl91kTSlJpjrW0TE280JNNrGQH3KaW\nKWKlqWbpWlyrtbXVnLxirWSqz0vts9If68g1noiic+yllk2yvArafEy1sY6qtHLP5UYI5qwu2gPg\n3QC+Ua/4WAjhFIB7AGyrH7YNwL1ZV3Qcp2mIlW0SAETkZgB/DeA5TC1CuAPA5wAcCiFcUj9GAAwV\nnwFgeHi4LLi/v3/ua+44DjZs2FBu9/T0TBt5y5HobQA2AfiTEMJTIvKXmJLjJSGEICLmX4rPfOYz\n0WdrBJxlJ0sTHXTwhS98odzetGlTtK/RH64lS5ZMk2os0bUcY4nO0syS6HpUm6VyIeMGBwfR29s7\nrX4FM5HoeqJEo+sCjSV6d3c3Tp8+nQykmQuJro/jZ13UvbW11Vz11VrhM1eis+zW+4r36pZbbsGz\nzz4bTSix3gkOdLEkul4FZseOHeX21772tWnnPPbYY/jIRz5ijqoX79yuXbtgkdPBDwI4GEJ4qv75\ne5jq4MdEZE0I4YiIrAEwYBViTfzXD5A7ECdP1B0rNzmBNeuMH5iOckslTNR/PLheuoPzHyudMCHl\nAtTlc/voTpdKiGHl8+ZzxsfHkwkTrRl6VuSW5eZslPO+6OA5s6x0Ha1kityp9R88blPtJuNnqMds\nUu4vy11pJaxIJcMEbHfjnLnJQghHARwQkevqX92FKbn+KIAt9e+2AHgk64qO4zSN3FH0PwHwUH0E\nfS+Af4+pPw7fEZH7AOwH8PH5qaLjOLMlq4OHEHYCuLXBrrtme2Er4QPvGxoaKrcvv/zy6DhrBUqm\nkMZjY2PmkjsajmjiOullgdietqSUnkzAEiyVuw2ITQArcisVeafL13I9tUyQ5bqybN/UqqxALDsL\n+7ajowPj4+PZbkSur34WbGZZK6WmovJEJMqvpm1wbh/e1u+iNjtT1049P/1crFz8Fh7J5jgVxju4\n41QY7+COU2EWLC96KhwVAI4fP15us42iXRbabcakZn/p/NUpnzsQ21GcMEDby1zHlStXRvtSrjwR\nia5nuYX4PJ7NBKRDbfXYQsq2tu7FGhfQdWQbkcc59DNKza7Tn61kGOzWGh4ejvaxDW4FcXEd+Rkt\nXrw4es/0O8H3k1rTDbBDSXOSbepEHBrPi+44zoVDVWcLh6o6jjP/NApV9V9wx6kw3sEdp8LMm0R3\nHGfh8V9wx6kw897BReQDIvKCiLxUz/zSFETkmyIyICK76LumppkSkT4R+ZmIPCciu0XkcwtUj04R\n+ZWI/KZejz9fiHpQfVpF5FkR+eFC1UNE9onI70Rkp4g8s4D1mNd0aPPawUWkFcBfAfgggBsAfEpE\nbpjPaxIPAviA+m4rptJMbQDwBNS89nlgAsCfhhBuALAZwH+s33+z6zEK4M4Qwh8AuBnAB0Rk8wLU\no+BzAPbQ54Wqxx+GEG4OIRTzLBaiHkU6tLdiKqHKnjmtRwhh3v4BuA3AP9LnLwH40nxeU13/SgC7\n6PMLANbUt9cAeKFZdalf8xEA713IemAqn96vAbx9IeoBYF39pb0TwA8X6rkA2AdgpfquqfUA0APg\nFdTHwuajHvMt0S8HcIA+H6x/t1BcFkI4Ut8+CuCyZl1YRK4EcAuApxaiHnVZvBNTiTkeD1MJPBai\nPf4HgP8CgEOxFqIeAcA/icgOEfkPC1SPqwAcB/A3dZPl6yKyZC7r8Xs7yBam/jw2xYUgIksB/D2A\nz4cQorQxzapHCKEWQrgZU7+gbxORjc2uh4h8GMBACGFH6pgmPpd31tvjg5gynd69APUo0qH9rxDC\nLQDOokE6tIupx3x38EMA+ujzuvp3C8Wxenop5KSZmgtEpB1TnfuhEML3F6oeBWEqI+7PMDU+0ex6\n3A7gj0RkH4C/A3CniPztAtQDIYRD9f8HAPwAwNsWoB6N0qFtmst6zHcHfxrABhG5qp4N5pOYSvW0\nUDQ1zVQ92+w3AOwJIfzFAtZjlYgUGXC7MDUO8Hyz6xFC+FIIYV0I4UpMvQs/DSH8cbPrISJLRGRZ\nsQ3gfQB2NbseoRnp0JowmHE3gBcBvAzgv8739ei62wEcATCOqb+U9wFYgakBnn4A/wSgd57r8E5M\nyavfAthZ/3f3AtTjJgDP1uuxC8B/q3/f1HqoOt2B1wfZmt0eVwP4Tf3f7uK9XIj2wJRX45n6s/m/\nAJbPZT08ks1xKszv7SCb4/w+4B3ccSqMd3DHqTDewR2nwngHd5wK4x3ccSqMd3DHqTDewR2nwvx/\nRbO63Ql6l9cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1d5578ae80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(n=(0, len(faces)-1))\n",
    "def display_faces(n=0):\n",
    "    plt.figure()\n",
    "    plt.imshow(faces[n].reshape((64, 64)), cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Mean and Covariance of a Dataset\n",
    "\n",
    "You will now need to implement functions to which compute the mean and covariance of a dataset.\n",
    "\n",
    "There are two ways to compute the mean and covariance. The naive way would be to iterate over the dataset\n",
    "to compute them. This would be implemented as a `for` loop in Python. However, computing them for large\n",
    "dataset would be slow. Alternatively, you can use the functions provided by numpy to compute them, these are much\n",
    "faster as numpy uses machine code to compute them. You will implment function which computes mean and covariane both\n",
    "in the naive way and in the fast way. Later we will compare the performance between these two approaches. If you need to find out which numpy routine to call, have a look at the documentation https://docs.scipy.org/doc/numpy/reference/.\n",
    "It is a good exercise to refer to the official documentation whenever you are not sure about something."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__When you implement the functions for your assignment, make sure you read\n",
    "the docstring which dimension of your inputs corresponds to the number of data points and which \n",
    "corresponds to the dimension of the dataset.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def mean_naive(X):\n",
    "    \"\"\"Compute the mean for a dataset by iterating over the dataset\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mean: (D, ) ndarray which is the mean of the dataset.\n",
    "    \"\"\"\n",
    "    N, D = X.shape\n",
    "    \n",
    "    mean = np.zeros(D)\n",
    "    for n in range(N):\n",
    "        for i in range(D):\n",
    "            mean[i] = mean[i] + X[n, i]\n",
    "    \n",
    "    mean = mean/N\n",
    "        \n",
    "    return mean\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def cov_naive(X):\n",
    "    \"\"\"Compute the covariance for a dataset\n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    covariance: (D, D) ndarray which is the covariance matrix of the dataset.\n",
    "    \n",
    "    \"\"\"\n",
    "    N, D = X.shape\n",
    "    covariance = np.zeros((D, D))\n",
    "    mat = np.zeros((N, D))\n",
    "    mean = mean_naive(X)\n",
    "    for i in range(N):\n",
    "        mat[i] = X[i,:] - mean\n",
    "    for i in range(D):\n",
    "        for j in range(D):\n",
    "            covariance[i, j] = covariance[i, j] + mat[:,i]@mat[:,j]\n",
    "    return covariance/N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.4001343 ,  0.43423554,  0.47628099, ...,  0.32141529,\n",
       "        0.31364669,  0.31045455])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean_naive(faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def mean(X):\n",
    "    \"\"\"Compute the mean for a dataset\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mean: (D, ) ndarray which is the mean of the dataset.\n",
    "    \"\"\"\n",
    "    mean = np.mean(X, axis = 0) # EDIT THIS\n",
    "    return mean\n",
    " \n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def cov(X):\n",
    "    \"\"\"Compute the covariance for a dataset\n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    covariance_matrix: (D, D) ndarray which is the covariance matrix of the dataset.\n",
    "    \n",
    "    \"\"\"\n",
    "    # It is possible to vectorize our code for computing the covariance, i.e. we do not need to explicitly\n",
    "    # iterate over the entire dataset as looping in Python tends to be slow\n",
    "    N, D = X.shape\n",
    "    covariance_matrix = np.cov(X, rowvar=False, bias=True) # EDIT THIS\n",
    "    return covariance_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.40013435,  0.43423545,  0.4762809 , ...,  0.32141536,\n",
       "        0.31364691,  0.31045464], dtype=float32)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean(faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.03256897,  0.03221722,  0.02860266, ..., -0.00627397,\n",
       "        -0.00599633, -0.00466273],\n",
       "       [ 0.03221722,  0.03582188,  0.03403026, ..., -0.00892454,\n",
       "        -0.00845002, -0.0067429 ],\n",
       "       [ 0.02860266,  0.03403026,  0.03782977, ..., -0.0121586 ,\n",
       "        -0.0113095 , -0.00946599],\n",
       "       ..., \n",
       "       [-0.00627397, -0.00892454, -0.0121586 , ...,  0.03519637,\n",
       "         0.03136033,  0.02796677],\n",
       "       [-0.00599633, -0.00845002, -0.0113095 , ...,  0.03136033,\n",
       "         0.03363069,  0.03122917],\n",
       "       [-0.00466273, -0.0067429 , -0.00946599, ...,  0.02796677,\n",
       "         0.03122917,  0.03254727]])"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov(faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.23923854,  0.07970797, -0.05669095,  0.09328563,  0.08262545],\n",
       "       [ 0.07970797,  1.09486781,  0.10168578,  0.08854109, -0.0128053 ],\n",
       "       [-0.05669095,  0.10168578,  1.17698659,  0.09787378,  0.017956  ],\n",
       "       [ 0.09328563,  0.08854109,  0.09787378,  1.00562478, -0.08340757],\n",
       "       [ 0.08262545, -0.0128053 ,  0.017956  , -0.08340757,  0.98028052]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.random.randn(100, 5)\n",
    "cov(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.23923854,  0.07970797, -0.05669095,  0.09328563,  0.08262545],\n",
       "       [ 0.07970797,  1.09486781,  0.10168578,  0.08854109, -0.0128053 ],\n",
       "       [-0.05669095,  0.10168578,  1.17698659,  0.09787378,  0.017956  ],\n",
       "       [ 0.09328563,  0.08854109,  0.09787378,  1.00562478, -0.08340757],\n",
       "       [ 0.08262545, -0.0128053 ,  0.017956  , -0.08340757,  0.98028052]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov_naive(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the `mean` function implemented, let's take a look at the _mean_ face of our dataset!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD3CAYAAAAwh5neAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnW2MXdWVpt/l+qCq7KLsMpg42MIMFG0hGJMWygAhre7Q\ntOhMBCiKECg9cUYoEyUziJEYDU5PNFL/GA1/0ppW0hpNCxIcNQ2DaDJYqKc7bjo9E0QnBBLAJjgU\nw4ewY1wGV5VdZdcne37UvZd1Vp311r7XVffCyXoky+fcfe4++3zsuu/aa+21JaWEIAiqybpONyAI\ngrUjOngQVJjo4EFQYaKDB0GFiQ4eBBUmOngQVJiz6uAicpOI/EpEXhORPavVqCAIVgdp1Q8uIl0A\nXgVwI4DDAH4G4I6U0i8BYHJyMhzsQdBGhoaGxH52Nr/gnwTwWkrp9ZTSHIBHANxyFvUFQbDKnE0H\nvxDA22r/cO2zIAg+JHS34ySjo6PtOE0Q/MYxMjJCy8+mgx8BsF3tb6t9Vsr9999f2D///PNLtwFg\neHi4dHvjxo2F4zZs2NDY7u/vL5T19vY2tnt6egAAv/71r7F9+/bCcV1dXaXbALBu3brSMv253Rcp\nmkF6v7794osvYteuXcuOzYGNmbz//vtNfe/gwYO44oorstuRe5w+lz2vbmN9+9ChQ9i5c2dpmd1e\nqWxxcbF0e2FhoXCc3q8fd+rUKQwODmJ+ft79ni6bnZ1tbM/MzBSOO336dGP75MmThbLx8fHG9vHj\nxxvbJ06cAAB8+ctfxoMPPtjYt9/Rx/7whz8E42wk+s8AjIjIxSLSC+B2APvOor4gCFaZln/BU0oL\nIvLvAPwdgC4A300pvbxqLQuC4Kw5Kxs8pfQ3AP4m59i6TG6cuLu7dNseq8uYNM4ts8cx6e2Vna1E\nr3+nFWnMJLpuhz3O+15XV5fbjmZMCK+NVkJ798OaR+zcun52vzXsmek2isiqvFeeeQcU32ltSto+\noftBM/en0I6so4Ig+EgSHTwIKkxb3GTAcomh961E9yQvkyW5Uo3Ryihx2f7Z1p8r81m7ctvEpJ89\nF5PG3vlyryX3ulYq8+6BHlEHirLcfscrs/vZMpmYfkyiM/MxJHoQBNHBg6DKRAcPggrTNhvc2hDa\n3mARZK3axV60k7XFvO/YdrTiCmPfm5ubyx5rsO41r4yR6wprZpxD49m+LApNR4ktLCycVRRa2b7e\n1hFotg7dptnZWbeNK53bI9eFxmxwO05l991zZx0VBMFHkujgQVBh2ibRWSQbkzAMLQWtlPLK7KSA\nXPmbazbkSvSpqamWou2Yu5G1wytbXFx0zYNcFxTgm0FMQtvnomW0t2332WSTXClfr6OnpwfT09OF\nMnud3rnZ+8fcfN47tm7dusJ+SPQgCJYRHTwIKkzHJDqL0vFGYa3M0nIpd0RzcnLSLWsl+m2lOjwJ\ndvz4cVd6W/lVNre9bJ/JfI0+1/z8vHvdTKLnSm8mr/V86vHx8cL+3Nxc6XG2jEWolcnwsmups2XL\nFoyPj9N59V6dzGSx8t2jmUg2zzSzxC94EFSY6OBBUGGigwdBhWmbDW7tSmZDeLaetcWYO6PMPhoc\nHMS7777rHmfxbCwWnWUps7OHh4dx+PDhgv2st88555xCHX19faXb9liWIMCz92dnZ10bnLmgWGQY\ns5+1m7Ket6yvrw/Hjh3DmTNnSo+zrk2938wsMU3Z2MiWLVtw7Nix7Bl1Xn0rnVvfO11mx1BWw5Uc\nv+BBUGGigwdBhemYRGf5w7Ts0nKPuRtyJh0MDg7i2LFjheO8SQe2jEVn6e8xd0b9HgwPD+Ptt98u\nyGsmw9evX9/YHhgYKJTpdNH6eyxyULfp9OnTrtxjEy30cwF8Sa1lNwBMT08Xzg0AO3bswDvvvFMo\n09+zdeS6yby8a0DR9ajvx4kTJ7ITUWhyJTPgmzpW1rNItnCTBUEQHTwIqkx08CCoMB1L+KCx9pG2\nsaytp2GuCM+tderUqcJxbNaSTUjgHcfcZGW5v3ft2oWjR48WbHBtS1s7W9vgehsoLt+ky6wdr21y\nbX/OzMy4dqa9Tu3ysu4vbT/re2zt56mpqWXbO3bswNGjRwvL/XguM4C/ExqWVES3Xx83NTVFZ+9p\nchNbsjBZb1xjbm6ucFzuLEJL/IIHQYVZsYOLyHdFZExEDqrPhkVkv4iM1v7ftLbNDIKgFXIk+oMA\nvgPg++qzPQCeSindJyJ7avv3skpYpI91x3gzcXIlkT2f/p6Vd2wlSW+GlK2DudD097RMPnXqlOta\nYpFb9tyeK8/KQi3Z9f2YnZ11XZZWojP3l5bl3jZQlPJark9MTLgSnd1vizcrT5slgG9WLSws0CQg\nOXUA/D56rjEbAZh7nYwVf8FTSv8XwAnz8S0A9ta29wK4NetsQRC0lVZt8AtSSkdr2+8AuGCV2hME\nwSoiObJXRHYAeDKldEVtfyKltFGVj6eUCnb45ORko+LR0dHVam8QBIqRkZHG9tDQ0DJ3SKtusmMi\nsjWldFREtgIYW+kL+/fvL+xrt5ANw/Psymbyp5fZ4Ndeey2eeeaZwnG5M9J0m1joZI699aUvfQnf\n//73C9et7XN9b4Ci+2twcLBQdu655za2h4eHG9vafQYUXW91e3RmZgZ9fX2Fduhrtravvm5tLwNL\nNnTZ9smTJwvHaZu8XscXvvAFPPbYYwXXlbb37dgIs5G9ZamtDV523M0334x9+/bRpa01rbw7QPHe\nleXvr78f7J2rZyb63ve+57YPaF2i7wOwu7a9G8ATLdYTBMEakuMmexjAPwH4LRE5LCJ3ArgPwI0i\nMgrg92v7QRB8yFhRoqeU7nCKbmjmRGzGGHMt5boscuW7jfDy2mT3tZRiLj9bh81Bruuzcr6OvRYm\nNb3kCjYaLjcRArvHbDkh3Q4tJ7VbDPBdfnNzc9n3m0Wa6fvDZtd5ySr7+voKdbJn7S2PtVKZxnvX\nWYIKIPKiB0GA6OBBUGk6NtmESbDcOlgecC86KzfwH/Dlrx1B1iPDbGKEnfSir1tLS+a6ZAkldB12\nJF5PbGHykZ1bH2uv0xtht/eKSXR9j71oL4BLdH1/9LmYhNZyfWFhgSbw8OR27nG2zVpqe/nZcvY9\n4hc8CCpMdPAgqDDRwYOgwnTMBs9d3lZv587eAYr2tLbtbF50/T1rL46Pjze29Zpm9jhvTS3Aj6w6\nfvx44ThtP9totaGhIXjo+6rt7mZmYDH3jIa5A701x1hyRpskQtvxbIlgOxtO47lVra2u79WmTR9E\nWR85cqQQOWhz1OtzlyXzWKkdZcfWYVGbrRK/4EFQYaKDB0GFaZtEb8b1kxu9xpIweHnB3nnnncJx\n+ntakgPAa6+9Vvo9Fq1m+fjHP97Y/tjHPtbYTim5UtNeM5uI4uWNY+aMdS96UX+5yx9ZvMlCtszm\neNOmj5cow5a9/vrrhTI9uUVfi5Xyuv0XXXQRAOArX/kKnnvuucY+sDwHnm6LNrnsc9H7ue5dXXdf\nX1/BnAk3WRAEy4gOHgQVpm0S3cpYPWpsy7xIIjZya+cdv/fee6VlR44cKRxnJZhG5wzT8p3lO7OT\nQfR87QsvvLCwrdulzQg7kcBb4qjs2DpsEkn9O7Ozs+ju7nblto3sY6uXetKeLaFkV1fVbdT3UZs5\n9loOHz5cKBsb+yA1gX4/bOSdvj/aSzExMYFt27Y19u1SV/rZ6zn31tOh74Gd+OPdb92m7u5u1xxg\ndVjiFzwIKkx08CCoMNHBg6DCdMwGZ0ujem4ymyCBLS3sLX9kI5O0XWztGh31pt0sNomBdols3ry5\nUHbZZZc1ti+55JLG9qWXXlqw73SknHWBaDt+48aNhTJ9PXqbzbKy0XXa1mNRbXqcgCWl0DZnbqIC\nG72nxxq0TQwUx03smIq2tfU4jB0b0XXo53LJJZfg4osvduvXNr4eo7HvlR4vss/Tm9mn6e3tLdw7\nO4ZgxzY84hc8CCpMdPAgqDBtk+hWmrEoIC2ntBS00Ugam1ZWSyQtcS+99NLCcVqiWxmnTQI9IcFK\ndC2XdOpiAIWoKC1dBwYGCnUyV5uW5Vai63N7ct3Wabe9ZY2sfNR1sjTEZSmaVyrbsmVLwe3kucyA\n4r26/PLL3TbqCT3WvNP3Xkv0yy+/HBdc8ME6Hkzaa9PPusJ0O6yc1mVefrxzzjmn8O5b92juCqvx\nCx4EFSY6eBBUmOjgQVBh2maDMztE28hA0SbXrgJrdzCb07qr6uzYsaOwr20q28adO3c2ts8777zG\ntnVZaHvRtkPblTYEV59bX7Mdk9AuJLskkXduW4dn4/f29hb2dRttaDDLz67rZ7PatJ2p69i8eXPB\njmU2pi7bsmWL2w49o9C6UfW59LjJtm3b6Huln4Vuh3WFsaQR+j3z7ndvby8N+bX33yN+wYOgwuQs\nXbRdRH4kIr8UkZdF5O7a58Misl9ERmv/b1qpriAI2kuORF8AcE9K6eciMgjgeRHZD+DLAJ5KKd0n\nInsA7AFwr1eJlRRaQloXg9735AxQlF3a7QEU3SL6uPPPP79wHIvc0jOEtMyy8jE3z7idceWZB/Z+\naNlppbdddqeOnSWn77+NZMvNY8ZWc/UkupXGXu7vgYGBwvm8fPIAjxLT5p41ZzSeubF+/fpCmXVP\n2fvvtYNFFXq57bTp19vbW7h3tv+smkRPKR1NKf28tn0KwCsALgRwC4C9tcP2Arg164xBELSNpmxw\nEdkB4BMAfgrggpTS0VrROwAucL4WBEGHEJYrrXCgyAYA/wfAf0kpPS4iEymljap8PKXU0MmTk5ON\nikdHR1exyUEQ1BkZGWlsDw0NLUsOmOUmE5EeAH8N4KGU0uO1j4+JyNaU0lER2QpgzK8BOHDgQGFf\n27fWftb2o7aBcpe2Bcpzc/f19RVcJ0DRRrThjHpf24T2OG/ZW6Boq9Ztu8suuwyvvvqqGz5q3SrM\nntP16+Os/ant+nrZG2+8gYsvvrhg/+s/+HasQWeg0TO1gGLGGz3Lio1J1O/byMgIRkdHXdvd3u9c\nN5zGjhmUZUvZtWsXXnzxxYINznKa6+PY7Dp7bi/xZN0Gv/LKK3HgwIGCTT4xMVGoo/4e33uvO+y1\n1H5autQ4AfAAgFdSSn+qivYB2F3b3g3giZXqCoKgveT8gn8KwL8CcEBEXqh99scA7gPwqIjcCeAt\nALetTRODIGiVFTt4SulpAF7i7xtyT8TyOjN548lYWweb+VSXeLOzs8ui5tgyO1oa6nOzZYDsdXpJ\nBrdu3eoua8Su00pQ3RbtwrF1eDnNu7q63PtopaW3VDFQvLZcF46uf3Bw0J1ZZV1t+rkwc8lzUQL+\n/RgYGKDJDr3nyUwn+8y8KD07m4zlhvcSZ1giki0IKkx08CCoMB3LycbKPPnUjEzRx9Zl20oSnS2z\nw6LV2HI/Xj7y8847r3SE3X4H4CPbnrS37fDyrtl7z8q8EWRbvzcqb+u0SzJ5I88sks3Kd+bR0HjL\nB23atIlep9dGlqecmXS6TJsbXV1drglnz82IX/AgqDDRwYOgwkQHD4IK0zYbnEWhMTst153Gltz1\n1pMCipFy1lZidrfXXmbT1tt4/PjxZWMBzKbSdjfLM86W981dbjb3Wmz92v73ZvKxOrq6ugrPwhsz\nAPzZWID/zHLfsU2bNtHlq70ZdWx9PdtGva+frb0fue1gxC94EFSY6OBBUGE6JtFZFJoXTdXM8kce\ndgI/c6t4kVUWJq+979kJJRobnaXrYGaKLst1fzXjvtSwiD32LLzoMmZ6WPdorjxlEr3suJmZGaxf\nv57W790TFmFoyzxTyrruWCReuMmCIIgOHgRVJjp4EFSYttng1r7NTd6Qa0fZ+jwbxdqH2gZirrDc\nzDf2OG+tr66uLtfdY11LbIZXrktxtWFhrLk2rB3jyHVLskQLbBaXpqysboPnjr2wddy8NjHYeXPf\nP0v8ggdBhYkOHgQVpm0SneXHZlFA2mXEIoKsRPIkKpNLTNbmSjV7HHP5ebm5mpHoHrnHWVo1RbyZ\ncdbl55lcKSXXTLFmlV0CymsHy+OusW7IXKnMZDkry2lHSin7vWLEL3gQVJjo4EFQYdom0XOTKQD+\nZAIrXXOlIIua07Qy2mn3WVSezSV25syZxr6W6Ox87Ny57bDb3veYWcXqz4kaK/tcn4+t3Om1CeCR\nYR7NmDPePbDvdyumTu4SWABffVUTv+BBUGGigwdBhYkOHgQVpmORbNpmYZP22eR+vZ/r4rJ2pYZF\niTE7mF2Ltru1nT09PV2wwdmSuKyN+rpz69BYW05fG1syiF0nm0mV+2x1fc0kU/BylbPZjDapJbPJ\nvXPnvsOsPuY2ZEtWM+IXPAgqTM7aZH0i8qyIvCgiL4vIn9Q+HxaR/SIyWvt/00p1BUHQXnIk+iyA\nz6SUpmqrjD4tIv8bwOcBPJVSuk9E9gDYA8Bd6tDKHub+8qRgM/mxtTxjUjA3wonV4bl37L7enp6e\nLuyzpA4sN5dXB4vs0+2dmZnJdv14q60CRfOD5Zr3coQvLCy40Y0sAQYzAZhE95JjLCwsZL8TzK3F\nVkfV9847bn5+Pjv/O2PFX/C0RH092J7avwTgFgB7a5/vBXBr1hmDIGgbkuOQF5EuAM8DuBTAn6eU\n7hWRiZTSxlq5ABiv7wPA5ORko+LR0dFVb3gQBEtrq9cZGhpaJj2yRtFTSosArhKRjQB+ICJXmPIk\nIvQvxQsvvFDYHxwc1A0rlOnUxnpbL2IPFPOa2bxdZfLsrbfewkUXXVQ4rhMSfX5+Hj09Pasi0b2l\nndi89/r22NgYtmzZ0pJE1x4AYMnkqNOsRN+yZQvGxsayr4Xla/PmiudI9KmpKWzYsKElid7MCqin\nT59ubE9NTTW2T548CQDYuXMnDh06VCg7ceJEoY7x8XEAwDe/+U23rUCTbrKU0oSI/AjATQCOicjW\nlNJREdkKYIx9l4XeMRuc5dguW3+s7Hz64bI83bkhqM08TGZj6X3m9mDJFD2b0yZ19BIhnjx5Mjvc\nlXVw/dLOzMy4dej69XOZnZ1d1lnr5Ca5sPutzKhjblQgP4yajVfofe/dse/HmtngInJ+7ZcbItIP\n4EYAhwDsA7C7dthuAE9knTEIgraR8wu+FcDemh2+DsCjKaUnReSfADwqIncCeAvAbWvYziAIWmDF\nDp5SegnAJ0o+fw/ADbknsjJFSw4m0Zl9y6LXvLJcaQPkRy2xaymTWd3d3VhYWHC/x1x5zFWoj9NL\nCdXPaY/r7+9v2HJ1mES37jWNlujaBmdLLWs5PT09XTAr2HNiyUI0+twsaYSug5l6dj/X1bsaEp3N\nlmREJFsQVJjo4EFQYT4Uk01yJTqTv2wFSi0Fc6WfbbOXuMHus1F0K9F1W7xJKbbNLKqL5Sorc5Nt\n374d7777LjxyV8UEipJdt5/lbtOS3I7K5y7b04qb036PeVnYxCLmBdH3xz5PL7rRbnsRb7YdjPgF\nD4IKEx08CCpMdPAgqDAds8HZbBvPpmXhl4xmbCyNbjOzh9i1eO61xcXF7FBYlpzRO7e93/re1e/H\n9u3bcfz4cfceMJdRjjsQ4MsCa2ZnZ90y+7munyW2YGG3mvp7ZMOHy/DcmfZ+eLa13dfPM2zwIAia\nIjp4EFSYji1dxCK3tDTxcmwB+fmy9feYBGNuISbDWY53zx2zbt06163FJpS0mmhBo+/Ne++911R0\nX1kdgH8PrLz2nouN3GKrvrKowtyVQW2SB2DpXWP56W1bWLQac3t6z8xK9NxluxjxCx4EFSY6eBBU\nmOjgQVBhOrY2GZvM7tkvbKnY3OR71h7KnT2VmwCP5S23IZbaPmVhpl57bVvYEsSeG256etrNY27b\nwZYx9kJL2TOzeDOk2LPNtcFZMo86GzZswMzMTHYbWRh1rpuMbee6khnxCx4EFSY6eBBUmLZJdDZh\nnbnQmETPjWKqy+Suri66VA9L5MBcM6xNbDlbL5EDSzJo8eRkrrkxPz/vJn9kUtWW6WsrSy5Rtm+3\nW5klxp4Fy6/myfdWly5iuc+ZvPai1RYWFtqTky0Igo8u0cGDoMJ0LJItN+GDhuULt5SNSg8MDNBV\nGlk72Iislqe5E2DY96wk19F8NrWwvh79PTZRhE1KyR3NX40kDPY7uSmPc0fzve/Y42wKb2am5Ep0\nJq9zEpqERA+CYEWigwdBhYkOHgQVpm02OMuLbu0XL/qLHcfcX9q2a8YGz849TWxVTa6dmrv+mD1W\nl7FINn1dPT09hbbkjiGwiL2y5BJlxzG3IYswZBGMXhmLZGPLaFm89yUn2eZKZTbRCXO1rdTOOvEL\nHgQVJruDi0iXiPxCRJ6s7Q+LyH4RGa39v2ntmhkEQSs0I9HvBvAKgHNr+3sAPJVSuk9E9tT27/W+\nbCWFjlCzstCTePa43Nxcuj4mw3Nlj4XJyTJ3Rk9PT7Zss/Xb69RLFOnvWSmv69D3gyVksDD3lJXb\nXnuZhNb3iuWhZ26s3IQPHmzyim0XS1qSW8a2WUKJVZXoIrINwL8EcL/6+BYAe2vbewHcmnXGIAja\nhuT8lRORxwD8VwCDAP5DSulzIjKRUqovKywAxuv7ADA5OdmoeHR0dNUbHgQBMDIy0tgeGhpaNnK7\nokQXkc8BGEspPS8iv1t2TEopiQj9S/Hoo48W9s8999zG9sDAQKGsv7+/sa0lqF3Uvq+vr/Q4oFyG\nDg8P48SJE4XjVkOiaznMUjvXt/v7+5ct1eNNsAGKaZOnp6cLZXpVT2+FT6B88sPNN9+Mffv2uRF7\nllyJziLvyp7Zddddh2eeeabwfPVxzIRbTYk+MjKy7McoV6Kz+eDWxNDPc2pqqrFdf7Zf/OIX8dBD\nD5WW2fofeOABckV5NvinANwsIp8F0AfgXBH5SwDHRGRrSumoiGwFMMYqYUnpmD3qJUUA4Npsdt/b\ntvs5SQHKYA/TvuDAUgfXnRHg637pfea6Y244lvxRd6DcmXK54ahsVpt1gXp12j+abE06DzZGk7t2\nnS1nMwzPdu09trx02b7HijZ4SukbKaVtKaUdAG4H8A8ppT8CsA/A7tphuwE8kXXGIAjaxtn4we8D\ncKOIjAL4/dp+EAQfIpqKZEsp/SOAf6xtvwfghtzv6uVlgaKdbWWs3s+VQWy2UK5Ez51xxEwFFllV\nZ/PmzZicnHRdRixPPJtVlBt51yq5cpjlTPNmss3OzrrjECwffm6Zleie6bewsNBSwgf2brbqJstd\nlpoRkWxBUGGigwdBhWnbZBMr0bUEsy4uLUeYy4XJIG/CympIdJZ/i0loXcfExES25G0lKip32R6b\naEHTzJJGnolh76m3WuzMzIw7Sp/rkrN16jJ7jd4SSouLiy2tWNrMKLr3vtgljXLzujHiFzwIKkx0\n8CCoMNHBg6DCdMwG15FcNgTVs7tZsgNmA+Xa4LlL4lr7hy0Z5Nlip0+fzrZxc5fLzV0y2OLZlexe\nWbyZfWxcQ9u+p0+fzp6Rpt8DFn3IEll4EXUsuYQtb9UG96IWm3GTsWWwNfELHgQVJjp4EFSYjkl0\nvW9nPmnJnpu7zZZ5Lhcrl3Lle+5KjyxZg5VgnhzOlXStou/N/Py8a4o0Y8547kBmGtjr8iR12YSd\nOmxSjecqBfhkpNxIttxlr9jkIe+4ubk5d4kjYHl/8ohf8CCoMNHBg6DCRAcPggrTNhvc2tnahrD2\nhM7mwRIElmVLqePZUdZWYjZ47mwyBrNNVyOBn97X94Ctn6Z5//333TY2E9bbyjpubNwhN9kGexas\nDq+N9rlYvPEWlpAhd/lgu+2Vle17xC94EFSY6OBBUGE+lBJd7zOJzqKdPKyUYlFXuckUcpckshLX\ni4TKnZFm62zm3Ppc7UyOoWEzvHIlOrsfufXlXjPgPyeW2ILJa90vrCT3ysr2PeIXPAgqTHTwIKgw\nbZPodiQxV6Ln5hxnS+SwdjBp70k8trJm7iqW3d3drjS0bdJlud4C+zkbKWcRal6drI0MNnrNUjtr\n9H4ryxPZ79nnkJtTLjf1N0v4oGW47RMxih4EASU6eBBUmOjgQVBh2maDW5tB29bMBtczy1jCh1xb\n2rajlXWumL1vy7TNrM/V39/vjiEwm5PNhmN496e3t5fanBpW5uV4Z+2zOcz1/dCJOO3Yi95nedE1\nbGwk130J+Ha3fa9yXVy6Ppt0sS1LFwVB8NEl6xdcRN4EcArAIoCFlNLVIjIM4H8C2AHgTQC3pZTG\n16aZQRC0QjMS/fdSSu+q/T0Ankop3Scie2r793pftpJCy6JWJTqbbKJhEp25oDxYkgErJ3Wb9bms\nRNffy5V0tk4W8Za7aiiLQstdoVNj74f3nb6+PjfJQzPP3Uv0wSLeNIuLi9nuL+bGYpFsXi70jqwu\nSrgFwN7a9l4At55FXUEQrAGSEyggIm8AmMSSRP8fKaW/EJGJlNLGWrkAGK/vA8Dk5GSjYruoehAE\nq8PIyEhje2hoaNmIZq5Evz6ldEREtgDYLyKHdGFKKYkI/Uvx1a9+tbA/ODjY2B4YGCiUbdzY+DuB\noaGh0u8AxRVK9RxyoDgKW5d41113HZ599tnCcWxE1oNJdFZWl4VXXnklDhw44Ervdkn0z3/+83j8\n8cfd49iEklxyJpvcfvvteOSRR9ZUojPvRp1Pf/rT+PGPf7wqEv3MmTON7ampqULZ9PR0Y3tiYqKx\nferUKQDAt7/9bdx1110FU1WnGQc+MGNfeumlZdehyergKaUjtf/HROQHAD4J4JiIbE0pHRWRrQDG\nWB0sQQBLKKc7rn3QubO49ItpZ7UxG9x7WVi4qC3Tf0B0O84555yWZsNZvHzkzbi7vHOz2V657jmG\nrqO3t9ft1Oy5MJclm03m/SGzSSjZsr25Nrh957w/2M3Y4Lk59Vd8q0RkvYgM1rcB/AGAgwD2Adhd\nO2w3gCeyzhgEQdvI+QW/AMAPan8ZuwH8VUrpb0XkZwAeFZE7AbwF4La1a2YQBK2wYgdPKb0OYFfJ\n5+8BuCH3RFZieMvI2mO9mTdAvpssN684cwvpNjJZyNxkGrtkMnNPManp5eZmYwFW4ubOBGOyvJWZ\nYPqZ9fe0ACfnAAAJLUlEQVT3u8k9mPmSG33I3GQ2icPZ5lMDiu+qldNeXnSWG87WkfvMIpItCCpM\ndPAgqDDRwYOgwrRtNhlLXsdsFL3djD/Us7+YDW7r8Oxn5iZjNqf9nK255cFs69wsH8ytl5vEMXcc\nwtbvhefaUFWWrSfXtmbjGtqWtkkoczOp5NrgzNXmJW5cXFx07fOyfY/4BQ+CChMdPAgqzIdColsX\nmjdZXs8ss3WwEE4m0W0iRK+MJX9ksOgvzw2Xm4OdlVlp7EnBnp4eOutKw5JeePKaRfbZ6/eSb7AE\nj7lLF7Eln6xbLDeRQ647jUl0b/kmu7x07lJUlvgFD4IKEx08CCpM2yQ6g01Eyc2PxUbRNVZyWbnK\n2lWHSfnc0U4muVo1AXLltSb3vllYbngWhebJ9+7ubtfEaGYWmye9WSSljU5jM/u8CDU22YTJd9Ze\nT8oDMYoeBAGigwdBpYkOHgQV5kNpg3vuDTajxtpYnnuKudNy7Rxrq7diK9m1uDQ52UfKztdKPvKe\nnh436SKz83KTXrBkDfY4L4EHawdzf7WSJWclGzw34QOzz72xAbZGWkSyBUGwjOjgQVBh2ibRmczM\nnUzA5JiVxt4kjNykhfbcGpZDjrm4bLSaF42Um6zB1rnaEj0nYWLZfm7ONLvtyfJmoiA995R97p68\nnp2dzZb2XuIGez7bfu+dsxKd5RwMiR4EQXTwIKgy0cGDoMK0zQZnM6nYsbm2GEuYyGad6TIWtqnb\nZNdS88I0LbqNc3Nz7vesfZWbUKKZENc61j2V6ybLfZ7sOHsu71mzxQdykyKy3OTWbZW79G/u2mEs\nh7nnprU2uCV72eiso4Ig+EgSHTwIKkzbJLqVj8yVoqUrkyIs0sdz8VjJpc/FTABNq3m6dcKKubm5\nlnKyMXLvKVuuqVWJnisZvediI7dyXWFMQmupbSW6J/Pn5uaylyRibcxd3jc3uQfLnceIX/AgqDBZ\nHVxENorIYyJySEReEZFrRWRYRPaLyGjt/01r3dggCJojV6L/GYC/TSl9QUR6AQwA+GMAT6WU7hOR\nPQD2ALjXPRGR6BZPorPRZRbJpsmd2FLW5hxYVJ6u//Tp0wXJruV6bvSXJTfxRG6EVCsJJAB/pJzV\nb3OQ5U7WsNLbk+VsVVl7HIt89NrYzOqfnrxuZrXZ3KWuc1YXHQLwOwAeAICU0lxKaQLALQD21g7b\nC+DWrDMGQdA2ZKWYVhG5CsBfAPgllhYhfB7A3QCOpJQ21o4RAOP1fQCYnJxsVDw6Orr6LQ+CACMj\nI43toaGhZfIuR4N2A/htAHellH4qIn+GJTneIKWURIT+pfja175W2GdypL+/v7Hd19fX2LZpk3WZ\nXa2zTMJ8/etfx3e+853CZ0yueqtd2nbo9g4MDLhl9e9dc801+MlPfrLqEp0Fh5RJy6uvvhrPPfec\nOye5mVxommYl+vXXX4+nn3561SX6mTNn3OPKJPo999yDb33rW3RVEq+NLODGogOlpqamGtv163r4\n4Ydxxx130FwF9WMPHjzongfI6+CHARxOKf20tv8Yljr4MRHZmlI6KiJbAYzRE2XOsgL4sjWaXNcM\ngy3Vo1+sVnOVe+c6c+ZM4UXVf6Ds/WDuEs+mzV0qitm+lly3Gcvn7bVxbm7OtWmZq4rN7MvtqCwf\nOftDyZYWYnj55XNn4TVzvhVt8JTSOwDeFpHfqn10A5bk+j4Au2uf7QbwRNYZgyBoG7nDxHcBeKg2\ngv46gH+NpT8Oj4rInQDeAnDb2jQxCIJWyergKaUXAFxdUnRD9olWIZKtmaQRHkzqsIQPzSQ/8LBu\nMm13awnZjET3YEkGrP3JJu1oWrXJNd65ZmZmXJs2d0JJvZ6yslyJPj8/n2X72jpyx0YAvlyW/ryZ\nlXA9IpItCCpMdPAgqDDRwYOgwnQs6WKuS6BVl5lno7AEBCzv+mq4yawNrtuo3WS5OccZbHad3p6e\nns72fbdig9vveC406zZkSR1YssPcpIjezLXFxUWarMG7j7m56y02Aabezs37z4hf8CCoMCuGqraK\nDlUNgmDtKQtVjV/wIKgw0cGDoMKsmUQPgqDzxC94EFSYNe/gInKTiPxKRF6rZX5pCyLyXREZE5GD\n6rO2ppkSke0i8iMR+aWIvCwid3eoHX0i8qyIvFhrx590oh2qPV0i8gsRebJT7RCRN0XkgIi8ICLP\ndbAda5oObU07uIh0AfhzAH8I4HIAd4jI5Wt5TsWDAG4yn+3BUpqpEQBPwcxrXwMWANyTUrocwDUA\n/m3t+tvdjlkAn0kp7QJwFYCbROSaDrSjzt0AXlH7nWrH76WUrkop1edZdKId9XRoO7GUUOWVVW1H\nfUWJtfgH4FoAf6f2vwHgG2t5TnP+HQAOqv1fAdha294K4FftakvtnE8AuLGT7cBSPr2fA/gXnWgH\ngG21l/YzAJ7s1HMB8CaA88xnbW0HgCEAb6A2FrYW7VhriX4hgLfV/uHaZ53igpTS0dr2OwAuaNeJ\nRWQHgE8A+Gkn2lGTxS9gKTHH/rSUwKMT9+O/AfiPAHR4VyfakQD8vYg8LyL/pkPtuBjAcQDfq5ks\n94vI+tVsx2/sIFta+vPYFheCiGwA8NcA/n1K6WQn2pFSWkwpXYWlX9BPisgV7W6HiHwOwFhK6XnS\nznY9l+tr9+MPsWQ6/U4H2lFPh/bfU0qfADCNknRoZ9OOte7gRwBsV/vbap91imO19FLISTO1GohI\nD5Y690Mppcc71Y46aSkj7o+wND7R7nZ8CsDNIvImgEcAfEZE/rID7UBK6Ujt/zEAPwDwyQ60oywd\n2m+vZjvWuoP/DMCIiFxcywZzO5ZSPXWKtqaZqmWbfQDAKymlP+1gO84XkXoG3H4sjQMcanc7Ukrf\nSCltSyntwNK78A8ppT9qdztEZL2IDNa3AfwBgIPtbkdqRzq0NgxmfBbAqwD+H4D/tNbnU+d9GMBR\nAPNY+kt5J4DNWBrgGQXw9wCG17gN12NJXr0E4IXav892oB3/HMAvau04COA/1z5vaztMm34XHwyy\ntft+/DMAL9b+vVx/LztxP7Dk1Xiu9mz+F4BNq9mOiGQLggrzGzvIFgS/CUQHD4IKEx08CCpMdPAg\nqDDRwYOgwkQHD4IKEx08CCpMdPAgqDD/HwK/qMTK8KRLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1d55207400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def mean_face(faces):\n",
    "    \"\"\"Compute the mean of the `faces`\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    faces: (N, 64 * 64) ndarray representing the faces dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mean_face: (64, 64) ndarray which is the mean of the faces.\n",
    "    \"\"\"\n",
    "    mean_face = mean(faces)\n",
    "    return mean_face\n",
    "\n",
    "plt.imshow(mean_face(faces).reshape((64, 64)), cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To put things into perspective, we can benchmark the two different implementation with the `%time` function\n",
    "in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.45 s, sys: 0 ns, total: 1.45 s\n",
      "Wall time: 1.85 s\n",
      "CPU times: user 4 ms, sys: 0 ns, total: 4 ms\n",
      "Wall time: 6.62 ms\n"
     ]
    }
   ],
   "source": [
    "# We have some huge data matrix, and we want to compute its mean\n",
    "X = np.random.randn(100000, 20)\n",
    "# Benchmarking time for computing mean\n",
    "%time mean_naive(X)\n",
    "%time mean(X)\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.35 s, sys: 112 ms, total: 3.46 s\n",
      "Wall time: 3.28 s\n",
      "CPU times: user 60 ms, sys: 0 ns, total: 60 ms\n",
      "Wall time: 31.8 ms\n"
     ]
    }
   ],
   "source": [
    "# Benchmarking time for computing covariance\n",
    "%time cov_naive(X)\n",
    "%time cov(X)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can also see how running time increases as we increase the size of our dataset.\n",
    "In the following cell, we run `mean`, `mean_naive` and `cov`, `cov_naive` for many times on different sizes of\n",
    "the dataset and collect their running time. If you are less familiar with Python, you may want to spend\n",
    "some time understanding what the code does. __Understanding how your code scales with the size of your dataset (or dimensionality of the dataset) is crucial__ when you want to apply your algorithm to larger dataset. This is really important when we propose alternative methods a more efficient algorithms to solve the same problem. We will use these techniques again later in this course to analyze the running time of our code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def time(f, repeat=100):\n",
    "    \"\"\"A helper function to time the execution of a function.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    f: a function which we want to time it.\n",
    "    repeat: the number of times we want to execute `f`\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    the mean and standard deviation of the execution.\n",
    "    \"\"\"\n",
    "    times = []\n",
    "    for _ in range(repeat):\n",
    "        start = timeit.default_timer()\n",
    "        f()\n",
    "        stop = timeit.default_timer()\n",
    "        times.append(stop-start)\n",
    "    return np.mean(times), np.std(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "fast_time = []\n",
    "slow_time = []\n",
    "\n",
    "for size in np.arange(100, 5000, step=100):\n",
    "    X = np.random.randn(size, 20)\n",
    "    f = lambda : mean(X)\n",
    "    mu, sigma = time(f)\n",
    "    fast_time.append((size, mu, sigma))\n",
    "    \n",
    "    f = lambda : mean_naive(X)\n",
    "    mu, sigma = time(f)\n",
    "    slow_time.append((size, mu, sigma))\n",
    "\n",
    "fast_time = np.array(fast_time)\n",
    "slow_time = np.array(slow_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAEGCAYAAADv6ntBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOX1+PHPmS0BglAFAQEBEfeKgggIuODGIu5WcSku\ntQWlbnXf935tbV35SatoQVFcChUrFlBEtCqCVZFFDSCrKKASCCSZ7fz+uHcmk5BlEmbJJOf9es1r\n7n2eu5z7ZJKT5869zxVVxRhjjMllnmwHYIwxxuwqS2bGGGNyniUzY4wxOc+SmTHGmJxnycwYY0zO\n82U7gFQoKiqySzKNMaaRa9WqlVRXZz0zY4wxOc+SmTHGmJzX5JJZYWFhtkNosKxtqmdtUzNrn+pZ\n21QvlW3T5JKZMcaYxseSmTHGmJzXKK5mNMaYVFJViouLiUaju7yt/Px8ioqKUhBV41O5bTweDwUF\nBYhUe9FitSyZGWNMJcXFxeTl5REIBHZ5W3l5eeTn56cgqsanctsEg0GKi4tp2bJlnbdlpxmNMaaS\naDSakkRm6iYQCNS7N5yxZCYiQ0TkaxFZLiI3V1F/gIh8JCJlInJ9XdY1xphsa/3celo/tz7bYTRZ\nGUlmIuIFxgFDgYOAkSJyUKXFfgKuAh6ux7rGGNOojB8/niOPPJLLL7+8Tutt2bKFZ555Jk1RNVyZ\n6pkdCSxX1ZWqGgSmAKclLqCqG1V1ARCq67rGGNPYTJgwgWnTpvH000/Xab2ioiImTJiQpqjqxvPt\n13i+/Toz+8rIXqAjsDZhfp1blu51jTEm51x77bWsWrWKc845h3HjxvHpp59y4oknMmjQIE466aT4\nzcbLli1j8ODBDBw4kKOOOooVK1Zwzz338O233zJw4EDuuOOOCttdvXo1ffr0YcyYMfTu3ZvLL7+c\nuXPncvLJJ9OrVy8+/fRTALZv386VV17J4MGDGTRoEG+++WZ8/aFDh3L00Udz9NFHM3/+fADef/99\nhg8fzq9//Wv69OnD5Zdfjmpmh8yVTOxQRM4Ghqjqb9z5i4C+qjq2imXvBopV9eFk100caNjutjfG\n7Kr8/Hzatm0bn2//0o8p3f73I/eodZkjjjiCmTNnsscee7Bt2zaaNWuGz+dj3rx5TJw4kQkTJnDr\nrbfSu3dvzjrrLILBIJFIhE2bNnHRRRfx3nvv7bTNNWvW0L9/f95++232339/hgwZwkEHHcQjjzzC\nzJkzmTJlCv/4xz948MEH2W+//Tj77LMpKipi6NChzJ49GxHB4/GQn5/PypUrGT16NLNmzeK///0v\nF198Me+99x7t27dnxIgR3HnnnRy3d3sAdnToknTbbNq0idLS0vh8jx494tM1DTScqUvz1wOdE+Y7\nuWUpXzfxwKtSWFhY6zJNlbVN9axtatbY2qeoqCitl9Mns20RIT8/n/z8fDZv3sw111zDypUrERFC\noRD5+fn079+fv/zlL2zatIkRI0bQvXt3tm7dGk84leXl5dGlSxcOP/xwAA488EAGDx5Ms2bNOOyw\nw3j44YfJz89n3rx5zJ49m/HjxwPOJfObN2+mffv23HDDDSxevBiPx8OKFSvIz88nEAjQu3dv9tln\nHwB69uzJ999/D24yq+54S0tLd6rbbbfd6Ny5c5XL1yRTyWwB0ENEuuEkovOA8zOwrjHG7LItl9T+\nzUbsSsbKy1b1B7uuHnjgAQYNGsTkyZNZvXo1p5xyCgDnnHNOvAd3zjnn8Mgjj9C1a9cat5WXlxef\n9ng88XkRIRKJAM5N45MmTdrpH5Q//vGP7LnnnnzwwQdEo1HatWtX5Xa9Xi/hcHiXjrmuMvKdmaqG\ngbHATGAZ8IqqLhGR0SIyGkBE2ovIOuA64HYRWSciu1W3bibiNsaYhmDr1q106NABgBdffDFevmrV\nKrp27cro0aMZNmwYS5YsoWXLlmzbtm2X9nf88cfz97//Pf691xdffBGPo127dng8HqZMmRJPfg1B\nxu4zU9UZqrqfqnZX1QfcsvGqOt6d/l5VO6nqbqra2p3eWt26xhjTVFx99dXce++9DBo0qEICmTZt\nGv3792fgwIEsW7aMkSNHsvvuu9OvXz/69++/0wUgybrhhhsIhUIMGDCAfv368eCDDwLwm9/8hpde\neokBAwZQWFhIixYtUnJ8qZCRC0DSrS5Pmm5s5/ZTydqmetY2NWts7VNUVESrVq3qtE46TzPmqthl\n+dFu+1dZX1Xb1NT2DeECEGOMadSS+V7NpI+NzWiMMSbnWTIzxhhTb5kc5aPGOLIdgDHGGLOrLJkZ\nY4zJeZbMjDEmBQpGHUvBqGOzHUaTZcnMGGMagWeffZaXXnop22FkjV2ab4wxjcCll16atm3Xdr9Y\nQ2A9M2OMaWBWr17NkUceyVVXXUW/fv0444wzKCkpAWDixIkcd9xxDBgwgIsuuogdO3YAzriJTzzx\nBN988w2DBw+usK2jjjoKgM8//5xhw4ZxzDHHcOaZZzqDAVcyZswYrrvuOk444QR69uzJ+++/zxV3\n38cRZ57LmDFj4svNmTPHeSzNyF/z6xtuobi4GICHHnqI4447jv79+3PVfQ/Gh8QaPnw4d911F4MH\nD6Z37958+OGHKW0z65kZY0wt6vJdWOVlC6pYpnji3Fq3s2LFCp555hkef/xxLr74YqZPn865557L\niBEjGDVqFAD3338/zz//PL/73e/i6+23334Eg8H4uI3Tpk3jjDPOIBQKceONN/Liiy/Spk0bpk6d\nyn333ce4ceN22veWLVuYPXs2M2bM4Pzzz2fWhPE8eec+HHPpaBYtWkTHjh3585//zL/+9S9ablzH\nI89NYty4cdx000389re/5aabbgLgdxecx3/mfcDJ+xwAQDgcZs6cOcyaNYuHHnqIl19+OblGTYIl\nM2OMaYC6dOnCoYceCsBhhx3GmjVrAFi6dCkPPPAARUVFFBcXc/zxx++07hlnnMG0adO49tprmTp1\nKs899xyFhYUsW7aM008/HWCnUe8TDRkyBBHh4IMPpm3bthzcY18ADjjgANasWcN3333H119/zckn\nn4wEywiGQvQZMBCAefPm8fjjj1NSUsLPmzZx4D77cLK73REjRux0PKliycwYY2qRTE8q1iOrvGx9\nx2as/EiV2GnGK664gsmTJ/PLX/6SyZMn88EHH+y07plnnsmoUaMYMWIEIkL37t1ZsmQJBxxwALNn\nz0563yKy0yNjIpEIXq+X4447jgkTJlT4Pq20tJTrr7+ed999l06dOvF/N11PaTC403bT8YgY+87M\nGGNySHFxMe3btycUCvHqq69WuUy3bt3wer386U9/4swzzwScBxdv3ryZTz75BIBQKMSyZcuAuo/i\n0adPH+bPn8/KlSsB2F5SwvLly+NPiN5jjz0oLi7m9bfn1Ps468p6ZsYYk0Nuu+02jj/+eNq0aUPv\n3r3jF15UduaZZ3LHHXfEn0UWCASYOHEiN910E1u3biUSiTBmzBgOPPDAOsfQpk0bxo0bx2WXXUaw\n2Hl22m333Mu+++7LqFGj6N+/P+3ataPXwXXfdn3ZI2BMnLVN9axtatbY2qc+j4BJ9WnGTKrt0vua\n6utbB/YIGGOMaXCS+V7NpI99Z2aMMSbnWTIzxhiT8yyZGWOMyXmWzIwxphKPx0Mw4f4okxnBYBCP\np35pyS4AMcaYSgoKCiguLo7fqLwrtm7dym677ZaCqNLH99UiAMK7t69zfZV1wTK8q75GNn9f43Yr\nt43H46GgoKoBwGpnycwYYyoREVq2bJmSbW3cuJHOnTunZFvpUjDFGZ+xeOg5da6vUBeN4nvv3+S9\n8jdkx/b4MtVtN5VtY8nMGGNySHX3s6VrvWR51q4k7x9/wbt8CQDhQ/viWzQ/LfuqiiUzY4wxu6zZ\nXZcjkQjRVrtTduFVRPocQ8HFx2Vs/5bMjDGmkUh376syz9qV5TPRKMHjTyd49m+gef2+99oVlsyM\nMcbUmXfp/8h//I74fMntTxLd9+CsxWPJzBhjTJ34PnqHvKf/iETKH+OSzUQGdp+ZMcaYZKninzGF\n/PH3IZEwwZOrvkoxGyyZGWOMSUpg8hPkvTwegLKRVxI8/8osR1QuY8lMRIaIyNcislxEbq6iXkTk\ncbd+kYj0Sqi7VkSWiMhiEXlJRBr28xSMMaYRCsyeivr8lF5xF6EhDadXBhlKZiLiBcYBQ4GDgJEi\nclClxYYCPdzXb4Gn3HU7AlcBR6jqIYAXOC8TcRtjmoaCUcfGrwQ0FcnPm+PT2rwFJTf8mXDfzF1y\nn6xMXQByJLBcVVcCiMgU4DRgacIypwGT1Hla6Mci0lpEOiTE2UxEQkBz4LsMxW2MMU1TOIz/7akE\npj0XLyq57UminbplMajqZSqZdQTWJsyvA/omsUxHVV0oIg8Da4ASYJaqzkpnsMYY05R5vvqCvOcf\nxbvu2wrlDTWRQQ5cmi8iv8DptXUDtgCvisiFqvpCVcsXFhbWus1klmmqrG2qZ21Ts4bePofffzkA\nn93+9M517nu6jiGV260p1l2t+/azhXR8+zV2X/wxAGWt27Lu5PPo/vITKd8nSdQB9OjRo8b6mEwl\ns/VA4miSndyyZJY5AfhWVTcBiMhU4CigymRW24EXFhYm3ThNjbVN9axtapZL7VNTnOk4huraZldH\n66jvcdRUd8jf7kJKtqN+P8HhFxAaPpJ2gTxwk1mq95nKz02mktkCoIeIdMNJUOcB51daZjow1v0+\nrS9QpKobRGQN0E9EmuOcZjweWJihuI0xpnGLRuKTUrKdcM9+lF14FbrnXlkMqu4yksxUNSwiY4GZ\nOFcjPquqS0RktFs/HpgBDAOWAzuAS9y6+SLyGvA/IAx8Bvw9E3EbY0yjpkrexEfjsyVX3Uek10AQ\nyWJQ9ZOx78xUdQZOwkosG58wrUCVd+Cp6l3AXWkN0BiT8zI90G66ZOQ4VAlMeQr/3DfiRZHeg9K3\nvzSzEUCMMaYJ8k9/nsB/XkG9Df46wKRYMjPGmF2Qizdc+2e9Rt7UZ1HxUDr69myHkxKWzIwxJk0K\nRh0bvyWgofDNm0He5CcBKLv0eiJHHpvdgFLEkpkxxjQhec8+DEDZBWMJHz0sy9GkjiUzY4xp7ELB\n+KRolLIzLyV00tlZDCj1Gsc3f8YYY3ZWVop/7hv433o5XhQc8itCp16UxaDSw5KZMcbkIN87rxPt\n2oNo5+4QyKtYuaMY/zuv45/5Kp5tWypUBc8bk5P3kdXGkpkxxuQI78L349P5kx4BQD0eoh32Jtpl\nv3hdiz+ci+zYDkCk2wEET72IZo/d5lQ2wkQGlsyMMaZWDeFmbM+aFeT//YH4fGjgyXhWFeL5bhXe\n9c4rRnZsJ7J/T4KnXkTk4N6NNoElsmRmjDEN3dYt5D92K1JWGi8qu/wWZyJYhmftSjyrvyF/otNb\n23Hb40T3OzQbkWaNXc1ojDENWThEsyfvwrP5ByL7HLhzfSCPaPcDCQ8+LV6UyURWPHFugxg+zJKZ\nMcY0YHkvPIH36y+Itt6D0qvuy3Y4DZadZjTGmAbK987r+N+djvr9lF51P/qLNtkOqU4y2WOznpkx\nxjRQeZMfB6DskhuIdq/iFKOJs2RmjMmKXBygN9MkEiE4bCThASdlO5QGz04zGmNMtoWCeFYX4i1c\njHf5knhxuGc/guf8JouB5Y6kkpmI5AF3AiOBPVS1lYicBOynqk+mM0BjjEm0y/d8qSJFP+FZu8K5\npH3tyniV56vPiR5w2K4HWZtoND7Z7L6xeFZ9jYRDOy1WOvp28HjTH08jkGzP7BGgI3AB8JZbtsQt\nt2RmjEmplN+kHA7HJ1uMPQ0p3lrlYvlP3c+OByZAQavU7LeKOHwfv43/zSnxIu/yxQBEOnYluu8h\nRHocTP4zDzmVzQvSE0cjlGwyOwPYV1W3i0gUQFXXi0jH9IVmjDG7KBrFt+A9Av+cEC+S4q1o8xZE\nO3Un0nkfop33If8ffwXAs2Uz+RP+7FwCn8pRM8pK8L/3Jv63XsHz08YKVSV/eIhI94OgRcvywlgy\nayAawn1ktUk2mQUrLysibYEfUx6RMcbsKlW8ixcSePVpvKu/qVC1/a+voLu3rZis3GSmzVrg+98H\n+Oa+Qfi4U1MWTovrzo33BqN7dSE47Lx47ytyaN+U7acpSzaZvQpMFJFrAUSkA/AoMKXGtYwxJgvy\nH7oO37LPAIi2bkPw9FHk/+MvAOgee1a7XtnF15H/1H3kvTiOyP490b26JL1P74L38Gxcj+eH9cjG\n7/D8sD5eJ8VbiXQ/kOApFxA57CjweBpc7yvXJZvMbgUeAr4EmgOFwNPAvWmKyxhj6s237DO0eQHB\nUy4gdMIZkJcPbjKrSbjf8YQWfYL/vzPJf+peSu58CvyBqheORvC/+VJ8ttmTd1W73R23PEp0/54N\ndsDfXDiNWJukkpmqBoFrgWvd04ubVVXTGpkxxtSBZ8Wy+HRw+PkEh4+s+D1UksouuhrvN1/iXbOC\nwGvPEBx5xU7LyObvyf/bg3i/WRQvCx8+gOieexFt1xHdsyPRPfeixY0XAGTmCskmLun7zESkObAv\nUAD0EPc/DFX9MD2hGWNyXcYenaJK3kvj4rPBX/22/ttq1pzSMbfT7P6xBP7zCpFD+lSo9n04m7xJ\njyIl24m22h1P0U8AlF7zQFVbMxmS7H1mv8a5BD8IlCRUKbB3GuIyxpikeRe+h7dwccq2F+1+EMEz\nLyXvtWfIe/rBeHneU/fh//gdAMK9BlJ66fUUjD09Zfs19Zdsz+xPwFmqOjudwRhjTJ0Fy8h7+W8p\n32xo+Eh8ixfg/eqLeJn/43fQQD5lF4wlfMzwBvsdWFOU7NiMQWBuGuMwxph68b89Dc+mDUQ6dk3t\nhj1eSn97G5pw43Kk2wHsuO8ZwseeknOJrKE8dyxdkk1mdwB/FZHcev6AMaZx27qFwPTnAaq8UGNX\n6R57Uvq7W+PzJbc/ibbvlPL9mF2XbDL7BjgV+EFEIu4rKiKRNMZmjDE1ypv2HFKynfAvjyTyyyPT\nso/IYUeVz/hsbPaGKtmfzPPAJOBlKl4AYowxWeFZ9y2+d99APR6C543Jdjgmy5Ltme0B3Kmqi1V1\nReIr2R2JyBAR+VpElovIzVXUi4g87tYvEpFeCXWtReQ1EflKRJaJSP9k92uMaZwCL49HNEr42BFE\nO3XLdjgmy5LtmT0HXITTO6szEfEC44ATgXXAAhGZrqpLExYbCvRwX32Bp9x3gMeA/6jq2SISwBmF\nxBjThPkWzUebtSB4xsXZDiWjarqIozFf4FGbZHtmRwLPuD2reYmvOqy/XFVXuqOJTAFOq7TMacAk\ndXwMtBaRDiLSCjgamADOaCSquiXJ/RpjGiDvJ3Pj0/lP3Il3wVwoK63zdoIjLkR3+0XqAjM5K9me\n2dPuq746AmsT5tdR3uuqaZmOQBjYBDwnIj2BT4GrVXX7LsRjjMkS37wZ5D37cPn8wnn4Fs5D8/IJ\n9xpI+Mjjdl4pHEKKt1Z4Dlm0bQdCJ56ZiZAzrin3sOor2bEZJ6Y7kBr4gF7A71V1vog8BtyMc7vA\nTgoLC2vdYDLLNFXWNtWztqlZVe1zeKW6tvPfptPslysss+7EX/GLJQto8d23+D96G/9Hb8fr/Fed\nha+kGG+wbKdtrx50KltWr6lxf6moS9d2a9tnU1Hb8ffo0SOp7VSbzETkIlV93p2+tLrlVPXZJPaz\nHuicMN/JLUtmGQXWqep8t/w1nGRWpdoOvLCwMOnGaWqsbapnbVOz2tqnx7774n99EnluIiu7YCx5\nk52H1Le+8AoU2L7xO3yfvItv/hy8a5xry/KKnEcmqseDttgNClri2eCcwGl76nm0rebG5RpjqWdd\nurbblD9Xqfy9qqlnNhLnknxwLv6oigLJJLMFOIMTd8NJUOcB51daZjowVkSm4JyCLFLVDQAislZE\n9lfVr4HjgaUYY1Im3QMCB6Y8ReA/r6DioezSGwgfPTSezGJ0z70InXIBoVMuiMez/U+T0YLdoFkL\n5xlgCbHm2ggcJr2qTWaqOixhuoqT2MlT1bCIjAVmAl7gWVVdIiKj3frxwAxgGLAc2AFckrCJ3wOT\n3SsZV1aqM8Y0cIH/vIJ6fZSOuZ1In2OTXk/bdUxfUKZRSXbU/M9U9fAqyheq6hHJbENVZ+AkrMSy\n8QnTClxZzbqfA0ntxxiTWQWjjuVwqujVRcLxSQ3kUfr7e4kcWvm6L2NSI9mrGfetXCDOA832SW04\nxpjGwj9neny65Po/E93/0CxG0/DYFYupVWMyE5HYTdKBhOmYrsCSdARljKmfjD0MszYlO/C/Xv4n\no6kmsuKJc52LHLIdSBNQW89sRTXTCvwXeDXlERljcp5/5qt4ttnYBiZzakxmqnoPgIh8rKozMxOS\nMSanbd1C4K0p2Y5iJzYMVOOW7E3TlsiMMUkJvPE8UlpCuGc/fF98nO1w0s4SYcOQ7NiMxhhTK9m0\nAf+c6agIwbN/k+1wTBNiT5ozponIxMUhganPIeEQof4nEN17p4ug0856SU2XJTNjTI2STYKetSvx\nfTQb9foInnVZ+gPLIEuSDV+yN01XNzZjGc7o9h+r6s4jgRpjmozAa08jqgQHn4q27ZDtcEwTk2zP\n7NdAf+AHnOTVCWgHLMS53wwROU1VF6YhRmNMDvB9/hGal09oxIXZDsU0QcleALIEuEFV91bVo1R1\nb+APwGc4ie0p4Ik0xWiMyRGhIeeirXbPdhimCUo2mZ0PPFmp7CngAndMxT8DB6UyMGNMbtGWrQgO\n/VW2wzBNVLLJ7AdgRKWy4cBGdzofCKUqKGNMjohG4pPBUy9yHtViTBYk+53ZVcCrIrIYWIvzEM1D\ngHPc+r7YaUZjmhz/m+UjfYSOO7VO69oVgiaVkh0BZJaI7IPzvLG9cB7l8qaq/hirB2alLUpjTPps\n3xaf9L37BuHjKp+EqZrvg5nkvfZ0eYE/kOrIjEla0veZuYnr+VoXNMbkDFm/imaP3R6fz//HXygL\nlhI6+Zwa1gLvFx+TN+GhtMVl4yiaukrqOzMR6SYiL4rIUhFZk/hKd4DGmHqIRmtdxLvwfZrfOwbP\nD+sqlOe9OA7/9Or/b/WsWEr+k3cj0SjB4efvcqjGpEKyF4C8CERxLse/qNLLGJNBBaOOjY/KUUFZ\nSXyy+Y0X4p8xBYqLdl4uGiUw9VmaPXEHUlpCqO/geFXppTegIuT9cwKB154B1QqryoY1NPvrzUiw\nlNDAIQTPubzGWIsnzrWelMmIZE8zHgwMUNXa/90zxmSFf9bU+LRn03fkvTyewNRnCfcdTOiE0+N1\n+Y/d5tzgLB6Cv/otoaHn4p8/B4DwMcMhkEfe3x8k8MYLUFYaX09+3kyzh29AircS7tmPskuuB5HM\nHaAxNUi2ZzYPODydgRhjylXb+6rO9m0EZrwUny255kHCvzwSCQXxf/Afmt89Ol7n+/wjtEVLSq//\nE6Fh5+2UkML9T6D0yntQr4/ArNfi5fkP34hn8w9Euh9E6ZV3gc+GdjUNR7KfxlXAf0RkGvB9YoWq\n3pnqoIwxdRN48yVkR3F8PnL4UUQOPwr5YR3+OdPxv/8W4l61GOm0D6VX34/uuVe124scMYjSqx8g\n/4k7kFAQAO+6lUQ77E3JdX+EvGbpPSBj6ijZnlkL4N+AH+ces8SXMSaLZMuP+Gf/s8o6bdeJ4Mgr\n2P5oeQ+r5M5xNSaymEjPvpT+ofyKxWjrNpTc8GcoaLXrQRuTYsneZ3ZJugMxxtSPf/rzSLCMcO9B\n+D59v+qFAnnl03XoVUUOLP92oeSmv6B7tKtvmMakVbXJTES6quoqd3qf6pZT1ZVpiMsYkwTZ+B3+\nuW+g4qHsrMuqT2YpoHt1Sdu2jdlVNfXMvgRautPLAQUqX7qkgDcNcRnT6KXiyc+Bac8hkQihgSej\nHbumJC5jclG1yUxVWyZMJ/vdmjEmQ5wnO7/tPNn59IuzHY4xWWVJypgcFfjnBESVkD3Z2ZjkLgAR\nkW7AA8BhQEFinfugTmNMhvk++y8asCc7GwPJ32f2IrACZzirHekLxxhTF6GTz876k52LJ86lsLCQ\nHlmNwjR1NpyVMWmUios8qqMtWhIcem7Kt2tMLsrYcFYiMkREvhaR5SJycxX1IiKPu/WLRKRXpXqv\niHwmIv/elTiMSbU6Dz1VX8EyfB+9Uz47fCS0aFnDCsY0HRkZzkpEvMA44ERgHbBARKar6tKExYYC\nPdxXX+Ap9z3mamAZsFuSMRvTKHjWrMA37038H86OD0kFEDrhzCxGZUzDkmwyqzycVV0dCSyP3WAt\nIlOA04DEZHYaMElVFfhYRFqLSAdV3SAinYDhOBehXFeP/RuTk5rdPRrvt1/F5yNd98O76htnJi8/\nZfuxx7SYXJep4aw6AmsT5tdRsddV3TIdgQ3Ao8CNlN/EbUzGFIw6lsPJzB98KfoJ/8xX4/Peb79C\nm7cg1P9EwscMJ9qlR2ZOaRqTY5J+hoOItAL2Z+dL8+ekOqhK+z0F2Kiqn4rIsbUtX1hYWOs2k1mm\nqbK22Vnsy+Lq2qam+mTr/EU/0u6jmezx+Qd4wqH4MqtOu4wtB/RC/QEIAoWFKdlfXSS7nn12qmdt\nU73a2qZHj+Suk032PrOLcb7zKqbipfkKVDtuY4L1VDw92cktS2aZs4BTRWQYkA/sJiIvqGqVN9fU\nduCFhYVJN05TY21Ts9rapqb6muoOmTcN339nIZEwAOHDB+D77L8AtDnzItrUYZux3mNVe6upLhk1\nHYN9dqpnbVO9VLZNsj2zB4CzVfWteu5nAdDDvfl6PXAecH6lZaYDY93v0/oCRaq6AbjFfeH2zK6v\nLpEZkyvk+/Iz6v55M1DxEOo7mNApFxDdu3uDOpVo36eZXJBsMvMBs+q7E1UNi8hYYCbOwMTPquoS\nERnt1o8HZgDDcAY13gHYY2dMoyM/bSLw+iR8896Ml4WOHkZw+Plo+05ZjMyY3JZsMnsIuF1E7qvv\njdOqOgMnYSWWjU+YVuDKWrYxF5hbn/0bk1XbtxF480X8s/6JhIKolN/iWXbZjVkMzJjGIdlkdi3Q\nHrhRRH4dR2jkAAAYLUlEQVRMrLCxGU1jkM6ROgBaXD8S2VEMQPiIoyk76zJa3DIqLfsypilKNpnZ\nd1TGJCMUxLO6EO+KpXhWLIsXy45iwgf1InjOb4nuc0AWAzSmcUr2PrP30h2IMemWlt5X8db4ZLN7\nx+BZvRxJuLQ+puTGh4kcfETq9muMqSDZS/Pvra4umeGsjGlMZMuPeP/3Ab6F7+P96rN4uXfFMlSE\nSMeuRLsfRKT7QeQ/9zCAJTJj0izZ04yVh7BqDxwDTEttOMbsmnR+9+Wf+Sq+hfPwFC5GVAFQT/mF\nHCU3Pkyk2wHQPGFcATeZGWPSq97DWYnIEGBkyiMypiEJlsUn814cB4D6/IQPOYLwEUcTPnwABVee\nCljvy5hsSno4qyrMAl5OVSDGNDSyaQP5T94Vnw/1OZbIEYMI9+wHzVpkLS67idmYnSX7nVnlIaua\n44zgsbaKxY3Jed5F88kf/wCyvfwCj7Kxd2cvIGNMjZLtmS3HGYdR3PkdwGeA3ShjMi6t94RFIwRe\nn4T/9UmIKuGe/fB98XHq92OMSalkk5lfVSNpjcSYBiD/rzfj+3IBKkLZWZcROuUCCi4ZnPE47FSi\nMXVTazJznxJdLCKtVbWstuWNyTnR8v/TfF8uQFu2onT0HUQOsQs6jMkVtSYzVY2IyDfAHsB36Q/J\nmMyQH3/AP28GvnnlD4OI7HMgpWPvQffYM+37t96XMamT7GnGycC/ReQxnCdAa6wi3Q/nNCbVvAvn\n4Z/7b7yLF8TvF4spufUx8AeyFJkxpr6STWZj3Pe7K5Un+3BOY7JHFc+a5fHZZk84g9aoz0+o9yDC\nx55Cs4eucyotkRmTk5K9abpbugMxJtVk8/f4PnoH34ez8X63Kl4e2asr4WOHEzrqRGjZOq0x2KlE\nYzJjV26aNqbBavbAVXi/WRSf14LdEHdQ4JIHnwOR6lY1xuQgS2Yma2q6Xyzpe8nCYTwrl+JbvBDv\nlwvixd5vFqH+AOFeAwgfdSKRQ/pQcNmJTmWKE5n1vozJPktmJmflP3Y73mWfISXbd6orvfxmwr0H\nZXXYKWNM5lgyMznFu3BefNr3vw8AiHboTPiQPkQO6UOzR24BIDxwSMr2WTxxLoWFhfRI2RaNMalm\nycykVcqGnlLF/+/J5L32TLyo9OI/EDnkCLRth13btjEm51kyMw1fsIy85x7G/+FsVCR+b1j4uBFZ\nDswY01B4al/EmOyRrT/T7KE/OIksL5/Sq+7PdkjGmAbIemZml6R1BHug2T2j8Wz+gejue1J67YNE\n9943LfsxxuQ2S2amYVGFbUXxWc/mH4h0P5DSq+5HW++Rtt3a5fXG5DZLZiY7QsH4ZOCFx/Fs2oBs\n/t55LystX6zf8ZRddiME8rIRpTEmR1gyM5mjiuebL/F/OBvfJ+/GiwOzp1ZcrHkLZIdz71jZ6Ntt\ntA5jTK0smZmMCEx9Ft+Hs/Fs2rBTXdl5Y4i27YC27UC0TXto0TL+XVwqE5mdSjSm8bJkZtIjFMQ3\nP6H39fokAKK7tyXc/wTCR51E89sucRYdem5WQjTGNB6WzExKSdFP+OZMx//u63iKfo6XhwaeTPio\nk4gceBh4vFmM0BjTGFkyMynhWV2If9Zr+D6eg4RDAEQ6d8e7dgUAZZffkrJ91XS60E4lGtM0Zeym\naREZIiJfi8hyEbm5inoRkcfd+kUi0sst7ywi74rIUhFZIiJXZypm4ygYdWz5d1jVaH7n5fg/mAmR\nMOFeAyi5+RFK7numxnWMMSZVMtIzExEvMA44EVgHLBCR6aq6NGGxoUAP99UXeMp9DwN/UNX/iUhL\n4FMRmV1pXbOLCkYdy+HUrWfjmzcjPq35zQkdPYzQCWeg7TqmPkBjjKlBpk4zHgksV9WVACIyBTgN\nSExIpwGTVFWBj0WktYh0UNUNwAYAVd0mIsuAjpXWNZkUCpL3whP4574RL9r+6Kv2uBVjTNZkKpl1\nBNYmzK/D6XXVtkxH3EQGICJdgcOB+ekI0tROftxI/pN34V25DPX7kZDz/ZglMmNMNuXMBSAiUgD8\nE7hGVbdWt1xhYWGt20pmmcbo8PsvB+Cz25/euc59r6ptYnUbZr9B16l/x7tjG2Wt9uDbs8dwwIT7\nq10v2e1WuW4sxgb0s2qqn5tkWftUz9qmerW1TY8eyT1JMFPJbD3QOWG+k1uW1DIi4sdJZJNVdSo1\nqO3ACwsLk26cxqqm46+pbt/JjyAaJXzwEYTG3E6nlq3BTWa1tWlV9bHv53Lhp2Gfm5pZ+1TP2qZ6\nqWybTCWzBUAPEemGk6DOA86vtMx0YKz7fVpfoEhVN4iIABOAZar61wzF27SpIhvW4P16Ed5vFsWL\nRaMER1xI8MxL7F4xY0yDkpFkpqphERkLzAS8wLOqukRERrv144EZwDBgObADuMRdfQBwEfCliHzu\nlt2qqjMwKZf/xJ14vl6EZ9uWnepKrrqPSO9Bddqe3fdljMmEjH1n5iafGZXKxidMK3BlFet9ANhI\ns2kgW37EN38Ovo/eiZf5Fs4DINrqF0T260l0/0PJe+FxgDonMmOMyZScuQDEpI7v/bfwffQ23qWf\nIRqtUFd66Q1E9u/p3CvmDvIbS2bGGNNQWTJrZKp98nPCM8Lyn3kIAPX6CB/an3D/E8j/f/cAED5m\neCbCNMaYlLJkloOqTVjVkA1ryH/y7vh85ICehPqfSPiIo6FgN6fQTWZ1Zd+JGWMaAktmjZzv43fI\ne+5hpLQkXlZyy2NZjMgYY1IvYwMNmwwLBcmb+Aj5T92HlJYQ6js42xEZY0zaWM+skWp2/1i8q75B\nfX7Kzh9LePCp+OfPqde27FSiMaahs2TWQNX1ezGCZXg2rInPeld9Q7RtB0qvvJtot/1TH6AxxjQg\nlsxyWOCVv+FZvxrPd6uQTRsQ1XhduNdASn9zE7RomdS2iifOdYaWSVewxhiTRpbMckUoiG/hPPzv\n/CteFHjzpfi0ejxE23eO985Kr7ovfp+YMcY0dpbMGjj58Qf8776B77038Wz9uUJd2ekXE+3YFd2r\nC9H2ncDnL38itCUyY0wTYsksi6r9XizhdGHzP4yMj9IR6dyd0PGnk/+PvwAQOuPiOu3PLuQwxjRW\nlswakmgU7+cfEpj+fHmZx0Ooz3GEjj+daI9DnB6Xm8yMMcY4LJk1BNEIvgXz8L/xAt61KypU7Xjk\nFbTV7lkKzBhjcoMlswag+a2XxC/ciP6iDaFhI8mb/ASAJTJjjEmCJbMs8az6pnx6wxqibdoRHH4+\n4UFDwR+IJ7O6su/FjDFNkSWzLPAu+4z8R2+Lz5dediPho04CX3I/DktYxhhTkSWzDPN++j75T92L\nhELxsvDRw7IYkTHG5D5LZhnkm/cWec/+GdEoocGn4Z/zerXLWu/LGGOSZ8kszWL3kpWdN4a8KU8B\nEDxtFMEzLq4xmRljjEmeJbNdlOyAwLFEVnbB7wmddFaaozLGmKbFklk6bd0Sn1SPh7LLbyF81Inx\nMjuVaIwxqWHJLEXkh3V41izHu3o5njXL8axejmfL5nh96dUPEDmsfxYjNMaYxsuS2a7YURyfbHHj\nhTtVa34zpLQEwBKZMcakkSWz+ohG8X3wHwKvPl1e1LoN0S77Et17XyLuu7bdi4JLBmcxUGOMaRos\nmdWRZ/kS8l54Au+3X1Uo3/HYa1mKyBhjjCWzJMSuWAwNPBn/BzMBpycWPG80+ePvr3Fdu8jDGGPS\nz5PtABo62fx9fNr/wUzU5yd4ygXseGgS4f4nZDEyY4wxMdYzq0o0gveL+fjfnY530fx4cfjwAZSN\nvAJt1zFeZj0vY4zJPktmCeTnzfjmzcA/9994ftoIgPr8EHbGUSy95oFshmeMMaYalsxcse/FYqJ7\n7kXouFMJDRoCLVtnJyhjjDFJyVgyE5EhwGOAF3hGVf+vUr249cOAHcDFqvq/ZNZNpfARRxMafCqR\nA3uBx75SNMaYXJCRZCYiXmAccCKwDlggItNVdWnCYkOBHu6rL/AU0DfJdVOm9Pf3pmOzxphGJKpK\nKArhqBJRCHiEPC84/5MbVaddIgp53sy0SaZ6ZkcCy1V1JYCITAFOAxIT0mnAJFVV4GMRaS0iHYCu\nSay7y4onzmXK8h0ULS3GI+ARECQ+DaAKUYUo6rzH550Pt1aYd8twLhn1esR5F/B4BK8406qgOK9o\nbNrdVkQhFFXC6vzShKMQUuddgIBX8HvA5xH8An6P4PcK6i4Tiiohd92QOx/V2LGBRyRhGrb87KfN\n1q2Ie8wecWL3iCDixqdKNNYW7ruq2z6o+55YBsGoEoxAWUTLp6PKjDWlAIzokk/A4xxLwCvxaa/7\nh0HcGCVhOqqUt0ul99i+46TiZOXjTjxG5+ehCT8L572oKMBuP/xcYduxSVVn+Uj8Vf6LHI06S3nc\n4D0JxxCLIXaMTnxSIWR121UrxRObjibsP/GzV51oLLYohN3PSViVpT+HATj4F774zyDgiX3GnGmR\nivGVxww/bw0QWP0j4agSdD9rwYgzHfs9iMUea9vCImefB7X24fUIPg/4xPk8e93PW+zzEoqq+/mB\nddsjALRr5sEjzudE3J+jV8p/vrH2lViZOx2Nx6EVYnI+U7H2cd7D0fKfZzihzapqYgHyvUIzn9DM\nK+T7nPlgMB//kh/i+4gkfF4St7PT57by9ivlhGi04t+iiMJPZVEAWvqlUnsrpU6zkecFDxV/B2Kf\nyfjnKuHvmvNzcNZt5hX3b1isvQWvx22bhL8z4UrHsuWSjqRbppJZR2Btwvw6nN5Xbct0THLdXdb6\nufWp3mQO8sO6bRnf6xurSzO+z7rzwcYd2Q4i7Za4Sa3ufPBj/X6OS7fUb58/lETrtV6q+D3OP5A7\nwuX/xJVElJJI5azkge31bdf62RaqPjOWRaDqdFy7nY+t+u14pOZ/rFKt0V0AUlhYWK9lFgyECWt8\n/BSSSj0Pif934/y3p/H/8Dwk/OdH+X/5kvBOwrZi/z0lTsf/S3fjiK0Lzn+pXom9Kz6P86WhV5yP\nUERxe17i9ErceXHXjb88Gp+O//dFwn/KFY63/LijJPzXSvkxQuKxa7zXFIufhHkAv0DAA36PEhD3\nj4A481GVeNzBaPnxhLT8F6FyvLhxxNrCGztGd76mMz2JPa6oe+y4xxrVij9XEo4t8Rgr9/Rw9xtr\nn9i01/28xI4hqlIxBsqPDyr2ZhMl9kgloaz8P2ut0MtEKrZ/5W3F28z9XMViFcHt/YvzH7b7cwm5\nP4+q4ou1pT/WqxKN/3xjn71Yu1TuXcc+i7FebFglYdpZJ/Y58Sd8bgIeTVhXKvx+xd4Te4BaoUwq\ntFcsjth0rG087rF4EsrK26z89yBRRKEs6rxKIxKfVhL+PsjO05V/PlWpKoV449up+PNP/BtSud2h\n/OeY2IuPTce3kbhuwt+NxB5lJKEHF2sbXzVtVNPf5dr+Zvfo0aPG+phMJbP1QOeE+U5uWTLL+JNY\nN662Ay8sLKx2mf9Lrs0arZrapqmztqmZtU/1rG2ql8q2ydTleguAHiLSTUQCwHnA9ErLTAd+LY5+\nQJGqbkhyXWOMMU1YRnpmqhoWkbHATJze8bOqukRERrv144EZOJflL8e5NP+SmtbNRNzGGGNyQ8a+\nM1PVGTgJK7FsfMK0Alcmu64xxhgTY3cFG2OMyXmWzIwxxuQ8S2bGGGNyniUzY4wxOU+0tjFUckBR\nUVHuH4QxxpgatWrVqtohEaxnZowxJudZMjPGGJPzGsVpRmOMMU2b9cyMMcbkvCaTzERkiIh8LSLL\nReTmbMeTKSLyrIhsFJHFCWW7i8hsESl033+RUHeL20Zfi8jJCeW9ReRLt+5xyfGnEIpIZxF5V0SW\nisgSEbnaLW/ybQMgIvki8omIfOG2zz1uubUPzgOHReQzEfm3O2/t4hKRVe5xfS4iC92y9LeP8yDI\nxv3CGdNxBbAPEAC+AA7KdlwZOvajgV7A4oSyPwE3u9M3Aw+50we5bZMHdHPbzOvWfQL0w3kaxFvA\n0Gwf2y62SweglzvdEvjGPf4m3zbuMQlQ4E77gfnuMVr7OMd0HfAi8G933tqlvG1WAW0qlaW9fZpK\nzyz+pGtVDQKxp1U3eqo6D/ipUvFpwER3eiJwekL5FFUtU9VvcQZ9PtJ94vduqvqxOp+ySQnr5CRV\n3aCq/3OntwHLcB4E2+TbBpyxUlW12J31uy/F2gcR6QQMB55JKG7y7VKLtLdPU0lm1T3Fuqlqp87j\ndQC+B9q50zU97XtdFeWNgoh0BQ7H6X1Y27jcU2mfAxuB2apq7eN4FLgR99muLmuXcgq8LSKfishv\n3bK0t0+je9K0qRtVVRFpspe0ikgB8E/gGlXdmnhavqm3japGgMNEpDUwTUQOqVTf5NpHRE4BNqrq\npyJybFXLNMV2qWSgqq4XkT2B2SLyVWJlutqnqfTMknnSdVPyg9uNx33f6JZX107r3enK5TlNRPw4\niWyyqk51i61tKlHVLcC7wBCsfQYAp4rIKpyvKwaLyAtYu8Sp6nr3fSMwDedrnrS3T1NJZva06oqm\nA6Pc6VHA6wnl54lInoh0A3oAn7inB7aKSD/3iqJfJ6yTk9zjmAAsU9W/JlQ1+bYBEJG2bo8MEWkG\nnAh8RRNvH1W9RVU7qWpXnL8jc1T1Qpp4u8SISAsRaRmbBk4CFpOJ9sn2lS+ZeuE8xfobnKtlbst2\nPBk87peADUAI57zzZcAewDtAIfA2sHvC8re5bfQ1CVcPAUe4H8oVwJO4N9zn6gsYiHNufxHwufsa\nZm0TP6ZDgc/c9lkM3OmWW/uUH9exlF/NaO3iHNM+OFcnfgEsif2tzUT72Aggxhhjcl5TOc1ojDGm\nEbNkZowxJudZMjPGGJPzLJkZY4zJeZbMjDHG5DxLZsbUQkRuFZFnal8ypfu8X0Q2i8j3SS6vIrJv\nuuMypqGyS/ONaWBEZG+ce266qDOKQjLrKNBDVZfXslxX4FvAr6rhXQw16/sxJsZ6ZsY0PHsDPyab\nyIwxlsyMiRORm0RkvYhscx8UeLxbfrc7/h4i8qSIFCe8wiJyt1u3l4j8U0Q2ici3InJVDftqJSKT\n3GVXi8jtIuIRkROA2cBe7vb/Uc36N4jIBhH5TkQurVQ3XJwHR24VkbWx+Fzz3Pct7vb7i0h3EZkj\nIj+6pzYnx4ayqqVdPCJys4iscNd9RUR2r24/Sf0QjKknS2bGACKyPzAW6KOqLYGTcR4yWIGqjlXV\nAlUtwBkS62fgdRHxAG/gDOPTETgeuEYSnpxbyRNAK5zhf47BGXvuElV9GxgKfOfu5+IqYh0CXI8z\nXmIP4IRKi2x3t9ca57lbY0Qk9iyoo9331u72P8J5+OEfgb2AA3EGfr07iXb5Pc4zpo5x1/0ZGFfD\nfoxJG0tmxjgiOE+7PUhE/Kq6SlVXVLewiLQF/gX8XlU/A/oAbVX1XlUNqupK4GmcwWgrr+t1y29R\n1W2qugr4C3BRkrH+CnhOVRer6nbcxBOjqnNV9UtVjarqIpzxOY+pbmOqulxVZ6vzgMRNwF8Tlq+p\nXUbjjL23TlXL3DjOFhF7tJTJOEtmxuD8QQeuwfmDvFFEpojIXlUtK86jY14DXlTVKW5xF5xTg1ti\nL+BWyh9CmKgNzpObVyeUrSb5hzPuRcUHGiZuBxHpKyLvuqcwi3CSTpvqNiYi7dzjXS8iW4EXYsvX\n0i5dcJ5zFjveZTjJr6pjNiatLJkZ41LVF1V1IM4faQUeqmbRJ4CtwO0JZWuBb1W1dcKrpaoOq2L9\nzThPMeiSULY3yT/PagMVnwG1d6X6F3EerdFZVVsB43FOJYJzXJU96Jb/UlV3Ay5MWL6mdlmLM8p5\n4jHnq/M8K7tM2mSUJTNjcL4bEpHBIpIHlAIlQLSK5X6HcwruAlVNrP8E2OZeLNFMRLwicoiI9Km8\nDXWe4PwK8ICItBSRLsB1OD2iZLwCXCwiB4lIc+CuSvUtgZ9UtVREjgTOT6jb5B7XPpWWLwaKRKQj\ncEPC8dbULuPdY+jiLttWRE6rYT/GpI0lM2McecD/4fSavgf2BG6pYrmROH+gv0u4ovFWN0GdAhyG\nc3/VZuAZnIs8qvJ7nAs1VgIf4PSmnk0mUFV9C3gUmAMsd98TXQHcKyLbgDtxkl9s3R3AA8B/3dOD\n/YB7gF5AEfAmMDVhWzW1y2M4PcBZ7r4+BvrWsB9j0sZumjbGGJPzrGdmjDEm51kyM8YYk/MsmRlj\njMl5lsyMMcbkPEtmxhhjcp4lM2OMMTnPkpkxxpicZ8nMGGNMzrNkZowxJuf9f6I5Lci9k7ZeAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1d48e8cac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(fast_time[:,0], fast_time[:,1], fast_time[:,2], label='fast mean', linewidth=2)\n",
    "ax.errorbar(slow_time[:,0], slow_time[:,1], slow_time[:,2], label='naive mean', linewidth=2)\n",
    "ax.set_xlabel('size of dataset')\n",
    "ax.set_ylabel('running time')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "## === FILL IN THIS, follow the approach we have above ===\n",
    "fast_time_cov = []\n",
    "slow_time_cov = []\n",
    "for size in np.arange(100, 5000, step=100):\n",
    "    X = np.random.randn(size, 20)\n",
    "    f = lambda : cov(X)               # EDIT THIS\n",
    "    mu, sigma = time(f) # EDIT THIS\n",
    "    fast_time_cov.append((size, mu, sigma))\n",
    "    \n",
    "    f = lambda :cov_naive(X)        # EDIT THIS\n",
    "    mu, sigma = time(f) # EDIT THIS\n",
    "    slow_time_cov.append((size, mu, sigma))\n",
    "\n",
    "fast_time_cov = np.array(fast_time_cov)\n",
    "slow_time_cov = np.array(slow_time_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAEGCAYAAADv6ntBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVOX1+PHPmbYVll5FUIOiUVGDwAKCFBtEESwRS1Bj\nYo3Rb2LvGn+WFKPRSKxBBVsUu0GN3YiiUhQRlipFpO+yy+5OO78/7t3Z2WXLADszO7vn/XrNa+7c\n59475z47O2eeW55HVBVjjDEmk3nSHYAxxhizuyyZGWOMyXiWzIwxxmQ8S2bGGGMyniUzY4wxGc+X\n7gCaQnFxsV2SaYwxLVxBQYHUV2YtM2OMMRnPkpkxxpiM1+qSWVFRUbpDaLasbupnddMwq5/6Wd3U\nrynrptUlM2OMMS2PJTNjjDEZr0VczVgfVaW0tJRoNBqbl52dTXFxcRqjar4yrW48Hg/5+fmI1HuB\nkzGmlWjRyay0tJSsrCwCgUBsXlZWFtnZ2WmMqvnKtLoJBoOUlpbSpk2bdIdijEmzFn2YMRqN1khk\npmUJBAI1Wt3GmNarRSezXdHu8TW0e3xNusMwxhizEyyZJdmUKVMYOHAgv/71r3dqva1bt/LII48k\nKaqd88Ybb3DPPfekOwxjjKlXiz5n1hw8+uijvPTSS/Ts2XOn1isuLubRRx/lvPPOS1JkiQmHw4wd\nO5axY8emNQ5jTPrkTz4SgNKp76c1joZYyyyJLr/8clasWMEpp5zCAw88wJdffslRRx3FEUccwdFH\nHx27YXDhwoWMGjWKYcOGMWTIEJYuXcott9zC8uXLGTZsGDfccMMO23766acZMmQIQ4cO5Te/+Q0A\nK1eu5Pjjj2fIkCGccMIJrFq1iuLiYg488MDYuaWysjJ++tOfEgqFmDp1KiNHjmTo0KGcddZZbN++\nHYALL7yQyy+/nNGjR3PjjTcybdo0rrjiCgDefPNNRo8ezRFHHMH48eNZv349AHfccQcXX3wx48aN\no3///kyZMqXBWDdu3MhZZ53FyJEjGTlyJLNmzUrSX8EY0xq0qpbZzpwLS2TZrec03Nq65557eOed\nd3j11Vfp2LEjJSUlvPnmm/h8Pt5//31uvfVWnnzySR577DEuuOACTj31VILBIJFIhJtuuomFCxfy\n8ccf77DdhQsX8uc//5m33nqLjh07smXLFgCuvPJKJk2axOmnn86TTz7JVVddxfTp0znooIP4+OOP\nGT58ODNnzmTUqFH4/X6OP/54Jk+eDMAf//hHpk+fziWXXALA2rVreeutt/B6vUybNi323oWFhbzz\nzjuICE888QT33nsvt99+O+Dczf/qq69SWlrKgAED+NWvfsWSJUvqjPXqq6/moosuorCwkFWrVnHS\nSSfx+eefJ/CXMcZkilS26FpVMku3kpISLrzwQpYtW4aIEAqFABg4cCB/+ctfWLt2Lccffzz77LNP\ng9v58MMPOfHEE+nYsSMA7du3B2D27Nk89dRTAJx22mncdNNNAEycOJEZM2YwfPhwXnjhhdihy2+/\n/Zbbb7+d4uJiSktLOfLII2PvMX78eLxe7w7vvWbNGs455xx+/PFHgsEgvXv3jpUdffTRZGVlkZWV\nRefOnVm/fn29sb7//vt89913sXW3bdtGaWkp+fn5iVeoMSbtmsshyJQlMxE5FrgX8AKPqOqdtcrP\nAK4CBNgGXKiq8xJZN1Fbz+lJRUVFg/dSVbXIGmt17Yrbb7+dI444gmnTprFy5Up+/vOfA3DKKacw\nYMAAZs6cySmnnMI999xDnz59mux9jzvuOG699Va2bNnCvHnzGD58OAAXXXQR06ZN46CDDmLatGl8\n8MEHsXXy8vLq3NaVV17JxRdfzNixY/noo4+4887qP0VWVlZs2uv1Eg6H640pGo3yzjvvZNR9bcaY\n5isl58xExAs8ABwHHABMEpEDai22HBihqgcBtwEP7cS6GaGkpITu3bsDMH369Nj8FStW0KdPHy64\n4ALGjh3LggULaNOmDdu2batzO8OHD+ell15i8+bNALFDdwMHDuSFF14A4LnnnqOwsBCA/Px8Djvs\nMK666iqOOeaYWIurtLSUbt26EQqFeP755xPehx49egDOubDG1BfrqFGjeOihh2LLzZ8/P6H3N8aY\nuqTqApCBwBJVXaaqQeAZYHz8Aqr6P1Xd4r6cBeyR6LqZ4ne/+x233norRxxxBJFIJDZ/xowZFBYW\nMmzYMBYuXMikSZPo0KEDgwcPprCwcIcLQPbff39+//vfM27cOIYOHcq1114LwN133820adMYMmQI\nzz77bI1W08SJE3nuueeYMGFCbN51113H6NGjOeaYY+jbt29C+3D11VczefJkRowYETt02JD6Yr3r\nrruYM2cOQ4YMYdCgQTz++OMJvb8xxtRFVJM/SLOInAwcq6rnua/PAgap6iX1LP8HoJ+qnpfIuvEj\nTccPKZCdnU3nzp13KtZuT28CYN2kxr+oTfpt2LCBioqKdIdhTIt26B+d+2TnXP9wk5UlKv6HdkMj\nTTe7C0BEZCTwK2DYrqwfv+PFxcU7nJNp7JxZMs6VZYrG6qY5atu2Lb169Ur6+xQVFSXcem2NrH7q\n15LqpqH92JWypqybVCWzNUD8N84e7rwaRORg4BHgOFXdtDPrGmOMab1Sdc5sNtBXRPYSkQBwGvBK\n/AIisifwInCWqi7emXWNMcbsuvzJR8Yusc9UKWmZqWpYRC4BZuJcXv+Yqi4QkQvc8inAjUBH4B/u\n+FRhVR1Q37qpiNsYY0xmSNk5M1V9A3ij1rwpcdPnAXV2RFjXusYYY0wV65uxlpbQ3DbGmNbGklkz\n89hjjyV0M3Jzc8opp7B169Z0h2GMaaWa3aX5rd25556b7hB2iqqiqgn3IGKMSY/m0odisljLLIlW\nrlzJwIEDufTSSxk8eDATJkygvLwcoN7hV+644w7+/ve/s3jxYkaNGlVjW0OGDAFg7ty5jB07lhEj\nRjBx4kTWrVu3w3uvX7+eM844g6FDhzJ06FA+++wzAO6//34KCwspLCzkH//4BwA333wzDz9cfVNj\nVQylpaWccMIJDB8+nCFDhvD666/HYhkwYADnn38+hYWFrF69moMOOohNm5y7KU4//XRGjBjB4MGD\n+de//hXbbs+ePbntttsYOnQoY8aMiQ0fU1+szz77bGxonMsuu6xGrynGGBOvVbXM8icfSaJ9sidy\n3iyRXzhLly7lkUce4b777uPss8/mlVde4Re/+MUOw688+eSTnH/++bH19t13X4LBYKzfxhkzZjBh\nwgRCoRBXXnkl06dPp1OnTrz44ovcdtttPPDAAzXe96qrrmLo0KFMmzaNSCRCaWkpc+fOZfr06bzz\nzjuoKmPGjGHo0KFMmDCBa665hrPOOguAl156iRdeeIHs7Gyeeuop2rZty6ZNmxgzZkxskM6lS5fy\n4IMPcvjhh++wzw888ADt27envLycUaNGccIJJ9ChQwfKysoYMGAAN9xwAzfeeCNTp07liiuuqDPW\nRYsW8eKLLzJz5kz8fj+///3vee6555g0aVKjdW6MaX1aVTJLh969e3PwwQcDcMghh/D9998DOw6/\nMnr06B3WnTBhAjNmzODyyy/nxRdf5PHHH6eoqIiFCxdy4oknAk7v8127dt1h3Q8//DA2QKbX66Wg\noIBPP/2UcePGxXrE//nPf86nn37KBRdcwMaNG1m3bh3btm2jXbt27LHHHoRCIW677TY++eQTPB4P\nP/zwQ6w11atXrzoTGcCUKVN47bXXAGfImKVLl9KhQwcCgQDHHntsrC7ee++9emN99tlnmTdvHiNH\njgSc3kk6deq0s9VvjGklWlUyK536fqNdNjX1ceXaw6JUHWasPfxKXYNwTpw4kcmTJ3P88ccjIuyz\nzz4sWLCAfv368fbbbzdJfFXGjx/Pa6+9xubNm2OdET/33HNs3LiRDz74AL/fz0EHHRTrB7G+IWI+\n+ugjPvjgA95++21yc3MZN25cbB2/3497D2GjQ8SoKpMmTYqNyWaMMQ2xc2ZpksjwK3vttRder5e7\n776biRMnAk4fZxs3boyNyhwKhVi4cOEO644YMYJHH30UgEgkQnFxMYWFhbz++uts376dsrIyXn/9\n9dgwMRMnTuSll17i5ZdfjrX6SkpK6NSpE36/nw8//JBVq1Y1ul8lJSUUFBSQm5vL4sWL+eKLLxpd\np65YR4wYwcsvv8yGDRsAZ+iYqlatMaZuzerWomBlSt/OklmaJDr8StXQLVUJJhAIMHXqVG666SaG\nDh3KEUccEUts8e68804++ugjhgwZwogRI1i0aBGHHHIIp59+OqNHj2bMmDGcddZZ9O/fH3CGaikt\nLaV79+5069YNgFNPPZW5c+cyZMgQnnnmGfbdd99G92vMmDFEIhEGDhzIzTffzIABAxpdp65Y+/Xr\nx/XXX8+ECRMYMmQIJ554Ij/++GOj2zKmpWtWCau2aATvgi/JevRu8n43MaVvnZIhYJItfgiYWvMp\nKCioMS/VhxkzSSb2ml/X3zgZWlLP58lg9VO/pq6bhr6jklHmf/NZsp55EIBohy5Ee/cl2vsnRHr3\nJdq7L3n/9wsAgsecgu+zd/Fs3bTDNur7Pt3ZusmoIWDSrTUmMWOMqYtn0XwCz/2z+vXm9Xg2r4c5\nn+ywbGCmc7ok2qUH4cIxhAaPJu+aySmL1ZKZMcaYHUjxZrL/cQsSjcbmld0xFe/KIjzuw7tyCVJW\nAkDwqJMIF44muvf+IPU2oJLGkpkxxpiaohGyHrwNz9ZNRPbrj3fRPAC0R2/CPXpD4RhnOVXyz3Zu\nnwme+dt0RQvYBSDGGGNqCbz4OL6Fc4gWtKfiohvrXzANLbD6tOhk5vF4CAaD6Q7DJEkwGMTjadEf\nYWNSzjvvMwKvPoWKh8oLb0TbdUx3SAlp0YcZ8/PzKS0tjd2oDM59UG3btk1jVM1XptWNx+MhPz/R\nDsqMMY2RjevI/uftAARPOpfI/oemOaLEtehkJiK0adOmxrz169fTq1evNEXUvFndGNO6ZT9wC1JW\nQrj/YELjTk93ODvFjtEYY4wBwLtsIdFOXan4zbWQYYfwMytaY4wxTUqKN8em1euj4uJbID9zTjdU\nadGHGY0xJhOloiciKd6M/41n8L/7cmxe5VmXEt27X9LeM5ksmRljTCsiWzfhf/NZ/O++jNTqDDg8\n8oQ0RbX77DCjMcakQbo6DM694nQC/3kOCVYSPmwo2295KOUxJIO1zIwxpoXzLJofm65KYsHxk4n2\naXwkjExhycwYY1qqaBT/608TePHR2KzttzzUopJYFUtmxhjTEpVsJfuh/4fv65rjHbbERAZ2zswY\nY1ocz6L55N5wHr6vP0fz21L+f3emO6Sks5aZMca0MDl3XoZEo0T6HkjFhTeiHbukO6Sks2RmjDEZ\nTLZsxLtwDt6Fc6rnRaMEx51OcOK54GsdX/OtYy+NMaaFyZp6D96Fc/D88P0OZeX/dxeR/oPSEFVN\nybzpuzY7Z2aMMRnCs3xRbNr/7st4fvgezc4h3H8wladdGCtrDoks1axlZowxGcA3679kPXJX7HXl\nSb8icsBhRPvsFzuUmPXMg+kKL+0smRljTHMWjRJ48TECrz5VY3bohLPSFFDzZMnMGGOaq/LtZD90\nO76vPkHFQ/D0i8ma9vd0R1VDKs+LNcSSmTHGNFM5f7wY7+rlaG4+FRffTOTAAWlJZs0lYTXEkpkx\nxjRT3tXLiXbfk/LL/h/abY90h9OsWTIzxpjmIBrB+9X/CLz5bGxWuP9gKi64HnLz0xhYZrBkZowx\n6VRZjv+j/+Cf+Tye9WtrFFVcdjt4vGkKLLNYMjPGmDTKu/wXSFkJANHO3QkdcwpZT93nFFoiS1jK\nbpoWkWNFZJGILBGRq+so7ycin4pIpYj8oVbZChH5WkTmisgXqYrZGGOSonx7bFLKSojscwDll9zM\n9rufInTUxDQGlrlS0jITES/wAHAUsBqYLSKvqOq3cYttBi4FTqxnMyNVdWNyIzXGmOTzrF0Zm95+\n/f1E+x6YxmhahlS1zAYCS1R1maoGgWeA8fELqOp6VZ0NhFIUkzHGpIVnQ/W5MUtkTSNV58x6Aqvi\nXq8GdqbzMAXeEZEI8E9VfagpgzPGmFSSDT+kO4QaMuE+ssZkygUgw1R1jYh0Ad4Wke9U9cO6Fiwq\nKmp0Y4ks01pZ3dTP6qZhVj/1q103vZZ8R1Y9ZQCHprgsnRqLp2/fvgltJ1XJbA3QK+71Hu68hKjq\nGvd5vYjMwDlsWWcya2zHi4qKEq6c1sbqpn5WNw2z+qlfXXWTXVkWm26o3lJdlmpN+blJVTKbDfQV\nkb1wkthpwOmJrCgieYBHVbe500cDtyYtUmOMSbL4c2ap0hIOJTYkJclMVcMicgkwE/ACj6nqAhG5\nwC2fIiLdgC+AtkBURC4DDgA6ATNEpCre6ar6n1TEbYwxTS4cQjZtSHcULU7Kzpmp6hvAG7XmTYmb\nXodz+LG2EqB/cqMzxpjUkE0/IhpNdxgtjo00bYwxKeRZ37yuZGwpLJkZY0wKSRrOl7UGCR1mFJEs\n4EZgEtBRVQtE5GhgX1W9P5kBGmNMS+JJ4j1mLf0ij4Yk2jK7BzgQOAPnBmaABcCFyQjKGGNaqto9\n45umkegFIBOAn6hqmYhEwbn3S0R6Ji80Y4xpecSSWVIkmsyCtZcVkc7ApiaPyBhjWirV3T7M2JoP\nJTYk0WT2PDBVRC4HEJHuwN9wOgw2xhiTiLISpLwMzc5FKrbXu5glrJ2X6Dmza4HlwNdAO6AIWIv1\nxGGMMQmruiw/2qV7miNpeRJKZqoaVNXLVTUf6Aq0cV9XJjc8Y4xpOaq6sdLOPdIcScuTcA8gIpIL\n/ATIx+lnEQBV/V9yQjPGmJZFqlpmna1l1tQSvc/sl8D9OBeClMcVKbBnEuIyxpgWp6plFu3Sw86L\nNbFEW2Z3Ayep6tvJDMYYY1qyqkE51VpmTS7RC0CCwPtJjMMYY1q8qhumo13snFlTSzSZ3QD8VUQ6\nJTMYY4xpscJhZNN6VATt1C3d0bQ4iSazxcAJwI8iEnEfURGJJDE2Y4xpMaqGftEOXcDnT3c4LU6i\n58yeBJ4AnqXmBSDGGGMSUNXzh13JmByJJrOOwI2qqo0uaYwxZgdVQ7+onS9LikQPMz4OnJXMQIwx\npiXz2D1mSZVoy2wgcImIXAf8GF+gqsObPCpjjGlhrPeP5Eo0mT3sPowxxuwCsX4ZkyqhZKaqU5Md\niDHGtGSx3j+sZZYU9SYzETlLVZ90p8+tbzlVfSwZgRljTItRtg3ZXopm50CbgnRH0yI11DKbhHNJ\nPtR/8YcClsyMMaYBnvVrALdV5nbSbppWvclMVcfGTY9MTTjGGNPyVF3JaH0yJk9Cl+aLyJx65n/R\ntOEYY0zLIxusT8ZkS/Q+s5/UniHOgGZ7N204xhjT8ljLLPkavJpRRJ5wJwNx01X6AAuSEZQxxrQk\nYlcyJl1jl+YvrWdagU+A55s8ImOMaWFi/TJ2tWSWLA0mM1W9BUBEZqnqzNSEZIwxLUgkjGz60Rn6\npWPXdEfTYiV0zswSmTHG7JpAyRYkGkXbdwZ/IN3htFiJXgBijDFmF2Rt2QCAWjdWSWXJzBhjkiiw\ndSNgF38kmyUzY4zZDfmTjyR/8pH1lmdtdVpmNvRLciXU0XADfTNWAquBWapa2WRRGWNMC5A/+Ujy\n3WkblDO5Eh0C5pdAIc5YZquBPYCuwBc495shIuNV1XoEMcaYOljLLLkSPcy4ALhCVfdU1SGquifw\ne2AOTmJ7EPh7kmI0xpiMZy2z5Eo0mZ0O3F9r3oPAGaqqwJ+AA5oyMGOMaUqNndtKJs3KRtu0S8t7\ntxaJJrMfgeNrzRsHrHens4FQUwVljDEtiQ39knyJnjO7FHheRL4BVgG9gAOBU9zyQdhhRmNMC1XV\noiud+v4urW/3mCVfoj2AvIXTQ/4UnPNk/wT2duejqm9VdX1VHxE5VkQWicgSEbm6jvJ+IvKpiFSK\nyB92Zl1jTOvR0OHCXT2UmOxDkHaPWfIl2jJDVTdRPfL0ThERL/AAcBTO1ZCzReQVVf02brHNOC3A\nE3dhXWOMabZs6JfkS3Rwzr1EZLqIfCsi38c/EnyfgcASVV2mqkHgGWB8/AKqul5VZ7PjubdG1zXG\ntCzpvFgjGaJdeqY7hBYv0ZbZdJwhYH4PbN+F9+mJc66tymqc82xNvm5RUVGjG0xkmdbK6qZ+VjcN\na8r6ObSBbaa6rCm2u7w8SKV9furU2Oemb9++CW0n0WT2U2CoqkYTXD5tGtvxoqKihCuntbG6qZ/V\nTcOSVT8NbTPVZY2u26snZOVUX7UYCcfK9jxsIASyGtx2a9SUn5tEk9mHOD8yvtzF91mDcwVklT3c\necle1xhjkiMSJvD8w7GX+eePRX1+NL8AzW8L2bnVy1oiS7pEk9kK4D8iMgNYF1+gqjcmsP5soK+I\n7IWTiE7DuRE7EbuzrjHGNL3SErL/cSu+BdU9+GkgCwlWIls3gttTvkmdRJNZHvAa4KdmKykhqhoW\nkUuAmYAXeExVF4jIBW75FBHphtPXY1sgKiKXAQeoakld6+5sDMYY0xQ83y8l+77r8Wz4gWibdni2\nbQWg7OGZUFmBlJYgpcVIaQk5d/8+zdG2HgklM1U9Z3ffSFXfAN6oNW9K3PQ6nEOICa1rjMlsu3sj\ncjp4Z79P9kN3IsEKIn32peLS28j7v19UL5CV7XRd1bFL+oJspepNZiLSR1VXuNN717ecqi5LQlzG\nmBYgExNWQ3LuvxmA0JCjqDznD3YurBlpqGX2NdDGnV4CKFC7czHFOfRnjDHNjyryw/f45n4am5Vz\n+6VEe/ZxHnvsRbRnnx3WYXspUrwZz5aNyNZN1UUeD8HTLiR09MnW12IzU28yU9U2cdM2IrUxJv0q\ny2OT3oVz0Lw2aH5bNK8tZGVXl309G++8T/HN/RTPhh9qbMK7eD7exfPr3HzuH05HijchwbrHGq74\nw5+I/PRnTbAjpqkl3J2VMcakmmzZiKfoG7xFX+Mt+gbPyuobbHPuvLzGsuoPVJf9+Yrq+W0KCB88\nGP8nMwEo/8Of8KxZgWfNcvexEqlw+oLwbFjrrJOVjbbrhLbvSLSgI/7P3gWwRNaMJZTM3MvibwcO\ngdgo4AC4A3UaY0zTiGt95V12co0i9VQfJIr06w+lJe7VgyVIKFhdtuc+RPoXEj6kkOje/cDjjSWz\nyEGHEzno8LiNKvlnjwSg7K4n0YKOkBN3jxjEkplpvlLVnZUxxjRKflxN9n3Vt65qdi6Rn/yUSN8D\nifY9kMg++5N//lgAyq+5t3pFVQhWkP+b45yy2x7diTetPvel3Xb6ziPTTLS47qyMMZnJO/d/ZP/z\ndmR7WWxe2YOvgieBa8xEnK6kTKuV6IUdVd1ZGWNM04pGCbz4ODn3XItsLyN82LDqskQSmTGkrjsr\nY0yGauxesfzJR3JoA+WNyb7nGnzzP0PFQ/CkXxEaN4n8c0bt0rZM65WS7qyMMaY+vvmfofltqbjw\nRiIHDkh3OE2qdOr7Ts/w6Q6kFUhZd1bGGENlBb6vPsb3yVuxWZE++1Lx21vRTt3SGNiuaym9m2S6\nhO8zE5ECYD92vDTfrlk1xjTIu+BLfP97C98XHyIV5TXKyq/7u3ULZXZboveZnQ08AJRS89J8Bert\nt9EY04rF3S8W33t8ZJ/9CQ85mqwn3UvrLZGZJpBoy+x24GRVfTOZwRhjWgZZv5bse6+PvY527k54\nyFGEhhwVu5crlsyMaQKJJjMf8FajSxljMlJT9m7v/fpzsh+8DSnbFpu3/U/TrWNek1SJJrO7gOtF\n5Da7cdqY9Gq2w6qo4n/jaQLPP4JolPAhQ/DN/Z9TluGJrNnVtdlBosnscqAbcKWIbIovsL4ZjWk+\n0pnosh64Bf9s532DJ04mOH5ys7pfrKE6sWSV+RJNZmcmNQpjTNLtbqILvPg40fad0Pad0A6dibbv\nBPkFsXL/7PfR7Fwqzr+WSHwvHk3MEo+pS6L3mX2Q7ECMMc2Pd87/YtOBl6fuUK5+f2w62r0X5Zf+\nEe3Re7ffd1cTliW61ivRS/Nvra/MurMypgVSxf/6dAL/fiQ2K3jCWciWje5jgzMK8/bSWPn2Gx+E\n3Py6tmZM0iV6mLF2F1bdgBHAjKYNxxgDab7II1hJ1qN345/135qzT/rVjstWlseGXbFEZtJpl7uz\nEpFjgUlNHpExrURzvCpRNm8g+77r8S5fhGbnUHH+deTE3S+2Axt2xTQTCXdnVYe3gGebKhBjTPrl\n3Hw+nuLNRDt3p+J3txPtZR38mMyQ6Dmz2p/oXOB0YFWTR2SMSa24m5s9xZsJ9zuEiktuhjbtkvq2\nzalFajJfoi2zJTj9MFbd+bgdmANMTkZQxpjkk00/4v/P8/g/eC02Lzj6RIKnXwK+3TloY0zqJfqJ\n9atqJKmRGGNSwvP9EvxvPIPvs3eRaM0OfYK/vCxNURmzexpNZiLiBUpFpJ2qVqYgJmNMHaR4M96v\nP4+9zv7r1eAPoP5A9bMrdkm9iNuVVHV3Urk3nAeAejyECscQOu4X5N7465TsgzHJ0mgyU9WIiCwG\nOgJrkx+SMaaKZ/HX+OZ/hnf+53hXLq5R5ps3q971Aq8+VW+ZBrIJHTmO0DGnJH1ATDsvZlIl0cOM\n04DXROReYDXO+TPABuc0pkmVl9VIUrm3/zY2rf4Akf0PxTf/M2fRy++AUCUSDEIoiISCZD11HwCV\nE891V1JQRVACLzk9eJTd82yNbqiMaQkSTWYXus8315pvg3Mas5ukZAverz7B9+VHeL/9CgmHYmXR\nrnsQPngQkYMHEul3CASyYvenRQ4p3GFbVcksNP6XO5RVJbOdTWTWujKZINGbpvdKdiDGtET13hgd\nrb6eKvfSiYg6BztUhMi+B+NdPB+A7XfXf7iwuSid+j5FRUX0TXcgplWz62+NSbVImKyH76x+7fUR\n/unPCB82jMhhQ9G27WNJ0BiTGEtmxuymneqWKhwi+8Hb8H3xYWxW2f0vQU5ecoIzppWwZGZMI/In\nH8mhNMG1M+6rAAAaRklEQVS5o2Al2fffhG/eLDQ3D9le5sy3RGbMbvOkOwBjmoP8yUcm99BeZQXZ\nf7vOSWR5bSm/6p7kvZcxrZC1zIxJgZy/XIV30TyibdtTceVfml0HvnbFosl0lsyMSQHvonlE23Wi\n/Kq/NMlIzMaYmiyZmVYjZeOHRcJ41qzAs+y72Kxox66UX/VXtGvP5L63Ma1UypKZO5jnvYAXeERV\n76xVLm75WJxe+c9W1a/cshXANiAChFV1QKriNqZBGusMh8C0v+NdtgjP90VIsGY3puXX3Yd27Jrw\nZpORcO1QomnJUpLM3M6KHwCOwukOa7aIvKKq38YtdhzQ130MAh50n6uMVNWNqYjXmHpFI3hWLcO7\naD7eRfPwLJofKwq89UL1Yl16ENmrH/7PnN7ediaRNaahpGQJy7RWqWqZDQSWqOoyABF5BhgPxCez\n8cATqqrALBFpJyLdVfWHFMVozI7Ky2KT2X+9Gu/ir5G4efEqJ55LdO/9iey1H+S3BYglM2NMcqUq\nmfWk5qjUq6nZ6qpvmZ7ADzh9QL4jIhHgn6r6UBJjNRlsd8+LeRZ/jXfFIjzLF+FdvghZV/2RrOoA\nONqpG5H9+hPZ72Ai+/Un76ozgbr7QzTGpEamXAAyTFXXiEgX4G0R+U5VP6xrwaKiokY3lsgyrVWm\n182h7nNd+7FDWSRM3ppltF32LfluWXwv9QBRjxdxx6VdceJ5lPbqS6igQ/UC28p37j3jXf9wVWGD\n+5QpMv2zk0xWN/VrrG769k2s189UJbM1QK+413u48xJaRlWrnteLyAycw5Z1JrPGdryoqCjhymlt\nmlPd7G4Lq6H9OGDF13i/+QLvd3OQivIaZZE99yHaZz8ie+1HdK9+RPfYi/zzjgag04Qz6bSL79lc\n6jVZmtNnp7mxuqlfU9ZNqpLZbKCviOyFk6BOA06vtcwrwCXu+bRBQLGq/iAieYBHVbe500cDt6Yo\nbtMSqOKdXz1Cc9UwKQDRHr0JHzggdvFG+W2Ppjw8Y8zuS0kyU9WwiFwCzMS5NP8xVV0gIhe45VOA\nN3Auy1+Cc2n+Oe7qXYEZzpX7+IDpqvqfVMRtMp934RwC/34U75JvYvNCg0YROfBwIgf+DO3QBah5\nJaIxJvOk7JyZqr6Bk7Di502Jm1bg4jrWWwb0T3qAJmMkcgjSs2QBgRcexfftVwBomwJkWzEAlRfd\nmOwQjTEplikXgBizU3Jvc34XaW4eweNOI3TUSeRfMDblcdh9X8akhiUz0yJ4ln1H4OWpsdeanUPo\n6JMJHnsq5LVJ6ntbwjIm/SyZmbSp73Chd/YHsemsf/2F8M+OILL/oeDz77ANz9JvCbw0Fd/8z2rM\nL/vT09C2XZPHbIxpniyZmeYjEibw/MME3nw2Nsv/3qv433sVzc0j3L+Q8M+Gxcqy/3QFvm9mA6BZ\n2YRGTyDwxtNOYRMmstKp7zuXEDfZFo0xTc2SmdktTdUTvRRvJvsft+D9bh7q9SIR50bl4Phf4v3y\nI7yrl+P/9B38n74TW8f3zWzncOLoCQSPOxXatKtOZsaYVsWSmUk7T9E3ZN9/M56tG4m260jFxTeR\ne/ulAAQnngsTz0V+XI3vy4/xffkR3iULnLLjzyR47CmQX7DbMdh5L2MymyUzk1SJtNxy7vgdEokQ\n2a8/FRfdiLbruMMy2nUPQmNPIzT2tNg2gyeft1OxWMIypuWyZGbSQjZvqJ6ORAgeeyrBU34DPvtI\nGmN2nn1zmEY15QjNsmYFgTefxfe/t2PzKi66ifCgkbu9bWNM62XJzKSEZ/F8Aq8/g2/u/wBQ8cTK\n6kpku5o47VCiMa2TJTMDOK2vQ9n1ZOB/ZwbqD4DP7zy701WqLuhQf4DwEccRPPYU8q48swkiN8YY\nS2atSpNdRr95A77P3q1xqDDryXsbXEfz2hAafSKhoyaibdvv1vtb68sYU5slM5Mw34dv4vv0bbwL\n5yCqNcpCo8ZDKAjhEBIKxqZ9C74EoOyvz0J2bjrCNsa0ApbMTJ1k83o8y77Du3RhbF72o3cBoD4/\n4UMKCRUeRc7fbwCgcvLldW6nqjVoicwYk0yWzFqrUBDZVoyUliClxbHZ2fdej2fZd3i2btxhlXC/\nQwgXjiF8+Iikd95rjDE7w5JZC1PfebH4+7ryzj8OqSivc33fVx8DztApkb32J7p3PwKvPgVAxTV/\na9JY7dyXMaapWDJr6cIh/DP/XWN4FKkoR71eNL8AbVOA5hfg+24uABW/uZbIPvujXXqCx7l8viqZ\nGWNMc2XJrAXzzv+MrGn341m3qsb80gdfg5w8EInNq2rRhYce3aQxWOvLGJMKnsYXMZko+97ryPnL\nVXjWrSLavRflf/hTdWFufo1EZowxmc5aZhmozvNiqni+XxJ76fvqEzQ7h+D4yYSOPqnOgS2bgrW8\njDHNgSWzTKaKZ81yfJ+9h+/z92scTgwNOYrgqeej7Tvt9ttYwjLGNHeWzDJY7rVn41m7MvY62qYd\nnm1bAag8/7odlm8oKdloysaYTGbJLIN4liwg8Nr06tdrV6L5bQn/bDjhQSOJ9OtP/rlj0hihMcak\nhyWzZip2Xuxf7+Gd/xmB15/Gu2hejWXK/3A3kf0PszHAjDGtnn0LplFCozDfcB7eVUsB0Nx8QqNP\njN33FTloYLJDNMaYjGDJLMkS7qk+GsWzdiWeoq/xLv4mNtu7ainRdp0IHXMyoZHHQ06e3cRsjDG1\nWDJLl/Ltscnsv16Nd8kCpGzbDotVnHsF4SFHOeODJcCuPDTGtEaWzFLEu3AOnhWL8axYjHflYmTd\n6liZb94sAKIdOhPpexDRvgeS9dR9AIRHjNthW5awjDGmJktmSSJrVuD/eGbsdc6dNYdIUa8PImEA\nKi64nsi+B6Edu8bKq5KZMcaYxlky2001zolt24p/1rv4Pv4P3hWLaywX2Xt/or37EumzL9E++xLt\n2Yf885x+EMOFdjm9McbsDktmTST73uvxzpuFuK0tzckjPHAk/g9eA6D8pgd3ant2KNEYYxJnyWxX\nRCN4F87FN+u/sVm+rz5GxUP44EGEhx1D+NChEMiKJbO6WMIyxpimYcksUap4ln6Lb9Z/8X3+Hp7i\nLTWKK39xAeEhR6HtOqYpQGOMab0smTWmbBv5Fx2/w+xo156EB48m8PITAITGnlbn6tb6MsaY5LNk\nVp+SrQRmPo//vy/FZkXbdyI8aBThwaOI9tkPRAhOPDeNQRpjjAFLZjuQzevxv/ks/vdfQ4KVNcq2\n//U58Nh4psYY09xYMnNVXWKvXl/sisTwIYUEjz+T6E9+msbIjDFNod3jawDYek7PNEeSeTKh7iyZ\n1RaNEBo0ktDPzyC65092e3PN6UPQnGJpbjKpbnY11qbYR1Vle1gpCyulIWVbKMqIVzYAuczvHqZd\nwEO+X/CI7FY8VWVbzu7BtpCysSLK+vII68uj/PK9zQA8M6YDPXK97JHnpX2WB4l7z6gqK7ZFmLcp\nyNyNIeZtDsXKjnh5PYO7BijsEmBw1yx65HkTjnXppG5srVS2BKNsrYyyNRjlvA+ci8FmHN2RTjle\nOmd76Jjtwe+pWQdRVUqCylZ33SNf3QDA3JO70i3HS7av8TrbHYnU9/LTu7OsJMzSkjDLqh7bwrHl\nhsz4kTbu37iN30Mbv/BkkdM1312DCmif5aFdwEO7LKF9locOWR46Zjdev01BVDU1byRyLHAv4AUe\nUdU7a5WLWz4W2A6crapfJbJucXFxwjtRVFRE3747DkFZ1TIru+tJtFuvhPerLqrKFxtCPL9sOw8t\nLAPgwA5+Ah7I8gp+j5DlgbfXOIcxz90vj3y/uA/ng3LpJ84gm9NHd0AAjwgeAY/ASW9tAuDlYzoh\nglvuPI99cyMAU0d2oCQYpSSkFAejlASjTPnWieV3B+bTKcdDlxwvXbI9dMrxcsTL6wGYNaEL67ZH\nWFsWYV15lB/KIjz8nbPeqXvn0CPPS888Lz1yneeeeV46Znvq/fJSVcIKFRGlMqJUhJXKiPN6qPue\nn0/oQpZXyPYK2T7nuesTawHnH09VUSCqzkOBgIcaX16N/T02VETZ95l1gPNPVxKMUhx06qY4GOWV\nlRUAnNA726kXt346Z3s4413ny3POSV3d91cnFkAVCl9y9mPuyV3xe4SAB/wewe8+V0ad/d4eVioi\nSnlYKY8ox73h/K0eHt6e8ohTXvX48zynn85z9sslHIWwQiTq1OWLy8tjf4+AV8jySuyz9bevSwG4\n6Kd5lIWcpLMtpJSFony8LghAYdcAbf1Cm4CHtu4X0r3fOOtd0b8NxUHnS7q40qmjWeud9dr4hdKQ\n87doiEegrV8oCHhYWRoBnM9/JKpEFMLufnzvlvXK9+IV8ABej+ABFhU7X6DZXqiINP43zvZCj1wv\ny7Y5C7f1CyWhxL4Weud7Y3Fef1hbNldG2FQRZXNFlE2VUTZVRGPlO6N9lrCl0omhXUAoDjZcd+0C\nQvdcL91yvby31vlu+HW/vNg6ivN5e2yR8/94Qu9strk/KEqCzvMP26OAUx8Bj+BzP48+j7C6zNmH\nfdp6iWj1/1NENbZestT3A6q+7+P6FBQU1PtPn5JkJiJeYDFwFLAamA1MUtVv45YZC/wWJ5kNAu5V\n1UGJrNsUyawpLNwS4t/LtvPvZeW79OHPZH6P848G7j+dOz+axI/X3m2cZLpHvs95zvNy2f+cHwG3\nDWjLouIwi7eGWVwcYmswNT/aWoMcr5AX9+PrG7fVs0eel+JglG0JJpFE5fqEztkeuuR46JTt5c1V\nzg+PMT2zWFsWYfX2CCV1/H275Xjo39HPwR0DHNLRH/tR8tpxnfh0XSWz1gf5fH1wp+JtFxDaZXmq\nWyABDzNWOD8uhnfPYmN5hA1uEqzrs9/GTfLtsqrrrWeulx/LI4SbwUc0zyfs3dbH3m297N3G5077\nGOf+SP54fBe2haJOq9z9W//O/Z87r18eW4NRtlQ6Lc8tldHYj4tUJLNUHWYcCCxR1WUAIvIMMB74\nNm6Z8cAT6mTXWSLSTkS6A30SWLdJfL6+kjL3l2fVr6Cq56j7a7zqF03Vr8yIwpqyCC8s3863W6qb\n491zPUzYK4d/LHB+RX14QmeCUQhGlFDUaZ2c+o7Twvrz4AJKQ85hmxL3gzJ9idN0P7ZX9g4tk3fd\nX21HdAvE5lfFWfUretye2RQEPLQNCG0DHgoCHq77vBiAG3/WlvXlETaUO4duNlRE+W6rE/s+bb10\nz635uMZd74Fh7VhbFmGt23JbXRZhgbvPoUZ+2OX7JNb6yvJCtlf41n3PvgU+KmKtNqf1Eqy1vfjW\nZ9U//bJtEfefJbjD+93wRUmN1/G/1M/rl0dBQNz68VAQEM553zlU9NiI9qyviLLBPaS1vjzCzNVO\nffdp48UDiDgt5aqYqupuz3wv4agTezCqhKOw3Q22U7aHbK+Q4xNy3OfP3L/VSXvlOPN9Qp77fMcc\np2X2l8ICfCL43F/XPoFfuYe1phzRnmDUqbNKt85u/dLZ79sGtCXf74klnTyfh/EznS+k147rFPsi\nKnGfb3HXu+qQNrRz66Sd+4Vd1dpfeUZ38nzOr/14VYeovjm1G+C0vKpavoe+8CPgfP59HsEr4BPB\n64FD/u2UzT25q/M/5rYSosDgGU5rd82Z3cnzN37R1baQcxRhoLved7/oRrfcmoe34r9Qh3XLApz/\n4wVbQgx/xTnkd9lB+XTM8tDBPVTYMcs58nCYux8rzuixw3s/Xkc8kahzKPInTztHA5ZO6kZBwLND\n3VWJqrKpIsq68ijrtkc45W3nu+GuQQUIzmcOnM//H2Y5/4//OrIDbQJCm7hDfm0DzucsrEowQuw5\nFFX6u/U9e2IXvHFHerwi/PQ5J87VZ3av84hHQ4emJ++XV29ZKqWqZXYycKyqnue+PgsYpKqXxC3z\nGnCnqn7svv4vcBVOMmtw3aZqmR3+4o8UFYfrLEtEu4Awvk8OJ++dy5CuAbz1fHDTpbFzJrvaag1F\nNdYqi/+nq9r7+uoh0XMmdf1zlYairCmLsMZNqqvLIqwujcR+BPx6/zz2LfCxb4Gf/dr56Jrjof2/\nqg9d7qxktuhbAquf+lnd1C8TDzOmLJkVFRXtcpy3Lg6wrlJi56Gqak3c4/keUfeZ2DF+j0COVxnW\nPkph+wgJ/Ig0xhiToPhk1xwOM64B4q+q2MOdl8gy/gTWjWksyzf0S+DJVv7jyX5B1s/qpmFWP/Wz\nuqlfU9ZNqtoRs4G+IrKXiASA04BXai3zCvBLcQwGilX1hwTXNcYY04qlpGWmqmERuQSYiXN5/WOq\nukBELnDLpwBv4FzJuATn0vxzGlo3FXEbY4zJDCm7aVpV38BJWPHzpsRNK3BxousaY4wxVexyBWOM\nMRnPkpkxxpiMZ8nMGGNMxrNkZowxJuOlrKPhZNqZHkCMMcZkpoZumraWmTHGmIxnycwYY0zGaxGH\nGY0xxrRu1jIzxhiT8VpNMhORY0VkkYgsEZGr0x1PqojIYyKyXkS+iZvXQUTeFpEi97l9XNk1bh0t\nEpFj4ub/TES+dsvuk0SHeW6mRKSXiLwnIt+KyAIR+Z07v9XXDYCIZIvI5yIyz62fW9z5Vj84Aw6L\nyBx3tA+rlzgissLdr7ki8oU7L/n1o6ot/oHTp+NSYG8gAMwDDkh3XCna9+HAYcA3cfPuBq52p68G\n7nKnD3DrJgvYy60zr1v2OTAYZ2ScN4Hj0r1vu1kv3YHD3Ok2OKOZH2B1E6sfAfLdaT/wmbuPVj/O\nPv0fMB14zX1t9VJdNyuATrXmJb1+WkvLLDbStaoGgarRqls8Vf0Q2Fxr9nhgqjs9FTgxbv4zqlqp\nqstxOn0e6I743VZVZ6nzKXsibp2MpKo/qOpX7vQ2YCHQE6sbwOkrVVVL3Zd+96FY/SAiewDjgEfi\nZrf6emlE0uuntSSznsCquNer3XmtVVd1htcBWAd0dafrq6ee7nTt+S2CiPQBDsVpfVjduNxDaXOB\n9cDbqmr14/gbcCUQjZtn9VJNgXdE5EsR+Y07L+n1k7Je803zpKoqIq32klYRyQdeAC5T1ZL4w/Kt\nvW5UNQIcIiLtgBkicmCt8lZXPyLyc2C9qn4pIkfWtUxrrJdahqnqGhHpArwtIt/FFyarflpLyyyR\nka5bkx/dZjzu83p3fn31tMadrj0/o4mIHyeRTVPVF93ZVje1qOpW4D3gWKx+hgIniMgKnNMVo0Tk\nKaxeYlR1jfu8HpiBc5on6fXTWpKZjVZd0yvAZHd6MvBy3PzTRCRLRPYC+gKfu4cHSkRksHtF0S/j\n1slI7n48CixU1b/GFbX6ugEQkc5uiwwRyQGOAr6jldePql6jqnuoah+c75F3VfVMWnm9VBGRPBFp\nUzUNHA18QyrqJ91XvqTqgTOK9WKcq2WuS3c8Kdzvp4EfgBDOcedfAR2B/wJFwDtAh7jlr3PraBFx\nVw8BA9wP5VLgftwb7jP1AQzDObY/H5jrPsZa3cT26WBgjls/3wA3uvOtfqr360iqr2a0enH2aW+c\nqxPnAQuqvmtTUT/WA4gxxpiM11oOMxpjjGnBLJkZY4zJeJbMjDHGZDxLZsYYYzKeJTNjjDEZz5KZ\nMY0QkWtF5JHGl2zS9/yjiGwUkXUJLq8i8pNkx2VMc2WX5hvTzIjInjj33PRWpxeFRNZRoK+qLmlk\nuT7AcsCvquHdDDXt72NMFWuZGdP87AlsSjSRGWMsmRkTIyJXicgaEdnmDhQ42p1/s9v/HiJyv4iU\nxj3CInKzW9ZDRF4QkQ0islxELm3gvQpE5Al32ZUicr2IeERkDPA20MPd/r/qWf8KEflBRNaKyLm1\nysaJM3BkiYisqorP9aH7vNXdfqGI7CMi74rIJvfQ5rSqrqwaqRePiFwtIkvddZ8TkQ71vU9CfwRj\ndpElM2MAEdkPuAQ4XFXbAMfgDDJYg6peoqr5qpqP0yXWFuBlEfEAr+J049MTGA1cJnEj59byd6AA\np/ufETh9z52jqu8AxwFr3fc5u45YjwX+gNNfYl9gTK1FytzttcMZd+tCEakaC2q4+9zO3f6nOIMf\n3gH0APbH6fj15gTq5bc4Y0yNcNfdAjzQwPsYkzSWzIxxRHBGuz1ARPyqukJVl9a3sIh0Bl4Cfquq\nc4DDgc6qequqBlV1GfAwTme0tdf1uvOvUdVtqroC+AtwVoKxngo8rqrfqGoZbuKpoqrvq+rXqhpV\n1fk4/XOOqG9jqrpEVd9WZ4DEDcBf45ZvqF4uwOl7b7WqVrpxnCwiNrSUSTlLZsbgfKEDl+F8Ia8X\nkWdEpEddy4ozdMy/gemq+ow7uzfOocGtVQ/gWqoHIYzXCWfk5pVx81aS+OCMPag5oGH8dhCRQSLy\nnnsIsxgn6XSqb2Mi0tXd3zUiUgI8VbV8I/XSG2ecs6r9XYiT/OraZ2OSypKZMS5Vna6qw3C+pBW4\nq55F/w6UANfHzVsFLFfVdnGPNqo6to71N+KMYtA7bt6eJD6e1Q/UHANqz1rl03GG1uilqgXAFJxD\nieDsV23/z51/kKq2Bc6MW76helmF08t5/D5nqzOelV0mbVLKkpkxOOeGRGSUiGQBFUA5EK1jufNx\nDsGdoarx5Z8D29yLJXJExCsiB4rI4bW3oc4Izs8Bt4tIGxHpDfwfTosoEc8BZ4vIASKSC9xUq7wN\nsFlVK0RkIHB6XNkGd7/2rrV8KVAsIj2BK+L2t6F6meLuQ2932c4iMr6B9zEmaSyZGePIAu7EaTWt\nA7oA19Sx3CScL+i1cVc0XusmqJ8Dh+DcX7UReATnIo+6/BbnQo1lwMc4ranHEglUVd8E/ga8Cyxx\nn+NdBNwqItuAG3GSX9W624HbgU/cw4ODgVuAw4Bi4HXgxbhtNVQv9+K0AN9y32sWMKiB9zEmaeym\naWOMMRnPWmbGGGMyniUzY4wxGc+SmTHGmIxnycwYY0zGs2RmjDEm41kyM8YYk/EsmRljjMl4lsyM\nMcZkPEtmxhhjMt7/B5DOi2ijfupZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1d51dff748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(fast_time_cov[:,0], fast_time_cov[:,1], fast_time_cov[:,2], label='fast covariance', linewidth=2)\n",
    "ax.errorbar(slow_time_cov[:,0], slow_time_cov[:,1], slow_time_cov[:,2], label='naive covariance', linewidth=2)\n",
    "ax.set_xlabel('size of dataset')\n",
    "ax.set_ylabel('running time')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Affine Transformation of Dataset\n",
    "In this week we are also going to verify a few properties about the mean and\n",
    "covariance of affine transformation of random variables.\n",
    "\n",
    "Consider a data matrix $\\boldsymbol{X}$ of size (N, D). We would like to know\n",
    "what is the covariance when we apply an affine transformation $\\boldsymbol{A}\\boldsymbol{x}_i + \\boldsymbol{b}$ with a matrix $\\boldsymbol A$ and a vector $\\boldsymbol b$ to each datapoint $\\boldsymbol{x}_i$ in $\\boldsymbol{X}$, i.e.\n",
    "we would like to know what happens to the mean and covariance for the new dataset if we apply affine transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def affine_mean(mean, A, b):\n",
    "    \"\"\"Compute the mean after affine transformation\n",
    "    Args:\n",
    "        mean: ndarray, the mean vector\n",
    "        A, b: affine transformation applied to x\n",
    "    Returns:\n",
    "        mean vector after affine transformation\n",
    "    \"\"\"\n",
    "    affine_m = np.zeros(mean.shape) # EDIT THIS\n",
    "    affine_m = A@mean + b\n",
    "    return affine_m\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def affine_covariance(S, A, b):\n",
    "    \"\"\"Compute the covariance matrix after affine transformation\n",
    "    Args:\n",
    "        S: ndarray, the covariance matrix\n",
    "        A, b: affine transformation applied to each element in X        \n",
    "    Returns:\n",
    "        covariance matrix after the transformation\n",
    "    \"\"\"\n",
    "    affine_cov = np.zeros(S.shape) # EDIT THIS\n",
    "    affine_cov = A@S@A.T\n",
    "    return affine_cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the two functions above are implemented, we can verify the correctness our implementation. Assuming that we have some matrix $\\boldsymbol A$ and vector $\\boldsymbol b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "random = np.random.RandomState(42)\n",
    "A = random.randn(4,4)\n",
    "b = random.randn(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can generate some random dataset $\\boldsymbol{X}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = random.randn(100, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that for some dataset $\\boldsymbol X$, the mean and covariance are $\\boldsymbol m$, $\\boldsymbol S$, and for the new dataset after affine transformation $ \\boldsymbol X'$, the mean and covariance are $\\boldsymbol m'$ and $\\boldsymbol S'$, then we would have the following identity:\n",
    "\n",
    "$$\\boldsymbol m' = \\text{affine_mean}(\\boldsymbol m, \\boldsymbol A, \\boldsymbol b)$$\n",
    "\n",
    "$$\\boldsymbol S' = \\text{affine_covariance}(\\boldsymbol S, \\boldsymbol A, \\boldsymbol b)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X1 = ((A @ (X.T)).T + b)  # applying affine transformation once\n",
    "X2 = ((A @ (X1.T)).T + b) # and again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One very useful way to compare whether arrays are equal/similar is use the helper functions\n",
    "in `numpy.testing`. the functions in `numpy.testing` will throw an `AssertionError` when the output does not satisfy the assertion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "correct\n"
     ]
    }
   ],
   "source": [
    "np.testing.assert_almost_equal(mean(X1), affine_mean(mean(X), A, b))\n",
    "np.testing.assert_almost_equal(cov(X1),  affine_covariance(cov(X), A, b))\n",
    "print('correct')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in the `???` below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "correct\n"
     ]
    }
   ],
   "source": [
    "np.testing.assert_almost_equal(mean(X2), affine_mean(mean(X1), A, b))\n",
    "np.testing.assert_almost_equal(cov(X2),  affine_covariance(cov(X1), A, b))\n",
    "print('correct')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check out the numpy [documentation](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.testing.html)\n",
    "for details.\n",
    "\n",
    "If you are interested in learning more about floating point arithmetic, here is a good [paper](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.22.6768)."
   ]
  }
 ],
 "metadata": {
  "coursera": {
   "course_slug": "mathematics-machine-learning-pca",
   "graded_item_id": "YoDq1",
   "launcher_item_id": "vCPZ0"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  },
  "widgets": {
   "state": {
    "12a53d600cc144b888c427ab187bb6df": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}